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Abstract 

Recent  progress  in  the  generation  and  sustainment  of  gas  discharges  at  atmospheric 
pressure  has  energized  research  in  the  field  of  plasma-aerodynamics.  Plasma  actuators 
are  promising  devices  that  achieve  flow  control  with  no  moving  parts,  do  not  alter  the 
airfoil  shape  and  place  no  parts  in  the  flow.  The  operation  of  a  plasma  actuator  is 
examined  using  a  macroscopic  (force  and  power  addition)  computational  fluid  dynamic 
model  of  a  dielectric  barrier  discharge,  DBD,  in  Fluent®.  A  parametric  approach  is 
adopted  to  survey  the  range  of  requisite  magnitudes  of  momentum  and  energy  delivered 
to  the  flow  field  and  to  identify  the  effects  of  this  localized  momentum  and  energy 
addition  on  the  flow  characteristics.  Simulations  consider  the  initiation  and  control  of 
flow  over  a  flat  plate  in  a  low  velocity  fluid.  The  simulation  velocity  profiles  are 
compared  with  the  experimental  observations  of  Corke  (AIAA  2002-0350)  as  well  as 
simulations  of  Font  (AIAA  2004-3574),  Boeuf  and  Pitchford  (JAP  97  103307  2005),  and 
Roy  and  Gaitonde  (AIAA  2005-4631).  The  simulation  is  extended  from  a  flat  plate 
simulation  to  examine  the  flow  modification  over  an  airfoil.  Flow  characteristics  of  lift 
and  drag  are  compared  with  experimental  results  of  Post  and  Corke  (AIAA  2003- 1024) 
and  the  compatible  energy/momentum  addition  is  identified.  Energy  and  momentum 
values  are  then  compared  and  related  to  characteristic  values  arising  in  DBD  operation. 
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MACROSCOPIC  COMPUTATIONAL  MODEL 
OF  DIELECTRIC  BARRIER  DISCHARGE  PLASMA  ACTUATORS 

I.  Introduction 

As  vehicles  are  pushed  further  and  further  along  the  envelope  of  powered  flight, 
certain  limits  are  being  reached  requiring  ingenuity  of  solution.  Specifically,  one  of  these 
limits  is  the  phenomena  of  stall  on  an  airfoil  or  lifting  body  when  it  is  flown  at  a  high 
angle  of  attack.  If  the  airfoil  is  forced  into  stall,  the  condition  where  lift  on  the  airfoil 
becomes  negligible,  the  vehicle  it  is  attached  to  has  a  tendency  to  either  fall  out  of  the  sky 
or  become  uncontrollable.  Generally,  the  stall  effect  occurs  when  flow  over  an  airfoil 
becomes  separated. 

In  the  past,  one  solution  that  was  explored  was  to  have  vacuums  either  sucking 
the  flow  back  to  the  airfoil  or  re-energizing  the  flow  by  blowing  into  it.  These  methods 
were  found  to  be  impractical  as  debris  eventually  clogged  the  tubes. 

Another  solution,  currently  in  use  today,  is  to  use  leading  edge  slats.  These 
devices  allow  flow  from  the  high  pressure  lower  side  of  the  wing  to  energize  the  flow  on 
the  upper  side  of  the  wing,  thus  preventing  separation.  However,  these  devices  cause 
unwanted  vibration  and  additional  drag  on  the  wing.  [9] 

More  recently,  experiments  have  proven  that  flow  can  be  reattached  and 
controlled  using  a  system  of  Dielectric  Barrier  Discharge  Plasma  Actuators,  which  will 
be  referred  to  as  DBD’s.  However,  the  mechanism  affecting  the  flow  is  not  fully 
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understood.  A  computational  model  of  the  system  is  needed  for  optimization  of  these 
devices. 

The  purpose  of  this  research  is  to  computationally  simulate,  evaluate,  and 
characterize  the  effects  of  the  addition  of  momentum  and  thermal  energy,  compatible 
with  the  operation  of  a  DBD,  to  the  neutral  gas  flow  over  a  flat  plate  and  an  airfoil. 

Background 

A  DBD  plasma  actuator  is  defined  as  “a  flow  control  device  with  no  moving 
parts,  does  not  change  airfoil  shape,  puts  no  parts  in  the  flow,  and  does  not  suck”  [1], 
The  basic  configuration  of  a  DBD  plasma  actuator  is  shown  in  Figure  1. 
However,  there  are  several  different  configurations  that  are  possible  as  seen  in  Figure  2. 

50pm 
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Figure  1 

Example  of  a  Dielectric  Barrier  Discharge  (DBD)  powered  by  an  AC  voltage  source.  [2] 
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Figure  2 

Varying  forms  of  Dielectric  Barrier  Discharge  (DBD)  configurations.  [3] 


Momentum  and  thermal  energy  additions  transfer  forces  to  the  flow  through  the 
combinations  of  collisions  of  electrons,  ionized  and  neutral  particles.  The  electric  field 
between  the  two  electrodes  causes  the  air  to  ionize  to  a  quasi- neutral  plasma  through 
acceleration  of  electrons  and  their  subsequent  ionization  collisions  with  air  molecules. 
With  each  collision,  given  the  electron  has  sufficient  energy  from  its  acceleration  by  the 
electric  field,  there  is  an  exchange  of  thermal  energy  and  a  high  probability  that  more 
electrons  will  be  freed  to  also  be  accelerated  causing  an  avalanche  effect.  This  allows  the 
flow  to  ionize  to  a  quasi- neutral  plasma  state.  The  same  electric  field  accelerates  the 
heavy  ions  in  the  opposite  direction  These  ions  transfer  their  momentum  to  the  neutral 
particles,  leading  to  the  modification  of  the  boundary  layer  flow  profile,  depicted  in 
Figure  3  and  Figure  4. 
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Figure  3 

Digital  Particle  Image  Velocimeter  (DPIV)  Flow  Velocity  Data 
of  a  DBD  in  Operation  at  -25  mm  [7] 


Figure  4 

DPIV  data  of  Boundary  Layer  Flow  Velocity, 

Normalized  by  Boundary  Layer  Thickness  and  Maximum  Velocity  [7] 
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Font  [8]  describes  a  single  breakdown  on  each  swing  of  an  AC  cycle  using  a  PIC 
code  employing  a  nitrogen  chemistry  model  Figure  5  displays  the  electrode 
configuration  of  the  modeled  system.  The  upper  or  top  electrode  is  exposed  to  the  flow, 
while  the  lower  or  buried  electrode  is  surrounded  by  a  dielectric  material.  The  simulated 
spans  of  the  electrodes  are  1  cm  deep,  while  the  buried  electrode  is  1.25  mm  wide  with 
the  exposed  electrode  0.25  mm  wide. 


The  results  obtained  suggest  that  the  majority  of  ionization  occurs  on  the 
“backstroke”  of  the  AC  cycle,  when  the  upper  electrode  goes  from  negative  to  positive. 
The  data  collected  is  displayed  in  Figure  6  and  Figure  7.  The  waveform  used  is  a  square 
wave  centered  about  ground,  0V.  The  first  part  of  the  AC  cycle  is  negative,  and  is 
referred  to  as  the  forward  stroke.  The  second  half  of  the  AC  cycle  is  positive,  and  is 
referred  to  as  the  back  stroke.  On  the  forward  stroke,  the  exposed  electrode  is  set  to 
-5000V  while  the  buried  electrode  is  kept  at  ground,  0V.  This  has  the  effect  of  causing  a 
few  random  electrons  to  start  a  breakdown  and  sending  the  free  electrons  onto  the  surface 
of  the  dielectric  barrier.  Within  30ns,  the  electrons  accumulate  on  this  surface  enough  to 
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nullify  the  field  between  the  two  electrodes  and  the  avalanche  ceases.  During  this  time, 


an  equal  number  of  ions  have  also  been  created  and  their  collisions  with  neutral  particles 
results  in  a  force  0.2  |iN  to  the  left  on  the  boundary  layer  flow  when  this  device  is 

operating  at  20W.  When  the  back  stroke  occurs,  the  exposed  electrode  is  set  to  +5000V 
while  the  buried  electrode  is  again  kept  at  ground,  0V.  The  electrons  on  the  surface  of 
the  dielectric  barrier  now  accelerate  towards  the  exposed  electrode,  again  causing  a 
breakdown.  However,  this  time  there  are  many  more  seed  electrons,  almost  all  from  the 
surface  of  the  dielectric.  Since  the  electrons  were  able  to  nullify  the  field  on  the  forward 
stroke,  there  was  -5000V  potential  at  that  dielectric  location,  which  results  in  a  total  back 
stroke  starting  potential  of  twice  the  forward  stroke  +5kV-  (-5kV)  =  +10kV  .  As  a  result 

of  more  seed  electrons  and  a  higher  starting  potential,  a  significantly  increased  amount  of 
ions  are  produced.  This  effect  can  be  seen  between 60- 70ns  in  Figure  7.  Because  the 
upper  electrode  is  exposed,  the  electrons  will  impinge  upon  it  and  do  not  nullify  the  field 
as  they  did  on  the  forward  stroke.  This  accounts  for  the  continual  increase  of  ions  in 
Figure  7  on  the  back  stroke.  The  ions  are  also  pushed  away  from  the  exposed  electrode 
and  result  in  a  force  1.4  p,N  to  the  right  on  the  boundary  layer  flow  when  operating  at 
20  W. 
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t  ■  61  MK 


t  ■  66  nuc 


Figure  6 

Charge  density  contours  during  Forward  (t=l-60ns)  and  Back  (t=61- 120ns)  Strokes  [8] 
From  left  to  right  and  top  to  bottom:  t=lns,  5ns,  61ns,  65ns,  10ns,  30ns,  70ns,  120ns 


time  (nsec) 


time  (nsec) 


Figure  7 

Computed  particles  during  forward  stroke  (left)  and  back  stroke  (right)  [8] 


Font  [8]  states  that  his  Particle -In- Cell  (PIC)  code  simulation  running  at  20  W 
with  a  voltage  between  1  to  5  kV  and  frequencies  of  1  to  10  kHz  will  produce  “...a  net 
force  of  6.0x10'  N”.  As  each  cell  is  6.25x10"  m  with  a  0.1mm  deep  span,  we  arrive  at  a 

unit  force  per  volume  of  — 6.0x10  N  _  qgQ^/m3  _  j  q5  jq/m3 

0.625xl0"V 
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Flows  are  more  likely  to  become  separated  if  they  are  laminar  than  if  they  are 
turbulent.  This  is  due  to  the  boundary  layer  being  much  larger  in  a  turbulent  case 
resulting  in  less  shear  force  in  the  boundary  layer.  This  effect  is  visualized  using  smoke 
flows  in  Figure  8  and  Figure  9.  The  flow  in  the  top  panels  of  both  the  curved  surface  and 
the  sharp  corner  surface  are  laminar  and  both  flows  are  seen  to  have  separated  boundary 
layers  near  their  highest  points.  In  the  bottom  panel  in  both  figures  the  flow  has  been 
“tripped”  to  turbulent,  thus  causing  the  boundary  layer  to  remain  attached  for  a  longer 
period  of  time  or  remain  attached  for  the  full  length  of  the  figure. 


Figure  9 

Flow  over  a  Sharp  Corner  Convex  Surface; 
Laminar  (top)  and  Turbulent  (bottom)  [6] 

The  DBD  operation  may  also  “trip”  the  flow  to  turbulent  earlier,  thus  maintaining 

attachment.  However,  due  to  the  low  velocities,  hence  low  Reynolds  numbers,  that  are 

under  examination  in  this  paper,  this  is  not  the  suspect  reason  for  maintaining  attachment 

in  reported  experiments  [7].  The  Reynolds  numbers  associated  with  the  velocities  under 

examination  are  <  0.5  x  106 ,  which  are  consistent  for  laminar  flows. 

By  adding  additional  momentum  and  thermal  energy  via  the  use  of  DBD’s,  the 

flow  is  energized  and  remains  attached.  Figure  10  and  Figure  11  depict  the  phenomenon 
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Figure  8 

Flow  over  a  Curved  Convex  Surface; 
Laminar  (top)  and  Turbulent  (bottom)  [6] 


of  energizing  the  flow  to  maintain  attachment  at  large  angles  of  attack  where  separation  is 
expected.  The  pictures  to  the  left  of  each  set  have  the  DBD  operation  set  off  and  show 
the  expected  separation.  The  pictures  to  the  right  of  each  set  have  the  DBD  operation  set 
on  and  show  attachment  being  maintained.  The  effect  of  energizing  and  maintaining 
attachment  of  the  flow  over  the  airfoil  will  increase  the  lift  coefficient  at  a  given  angle  of 
attach  as  well  as  increasing  the  stall  angle  for  the  airfoil. 


Figure  10 

Reattachment  of  Separated  Flow  with  Actuator  ON  for 
NACA  663-OI8  Airfoil  at  a=-16°  (Smoke  used  for  visualization)  [4] 


Plasma  off 


Plasma  on 


Figure  11 

Reattachment  of  Separated  Flow  with  Actuator  ON  for 
NACA  0015  Airfoil  at  a=12°  (Smoke  used  for  visualization)  [5] 
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Approach 

The  commercial  code,  Fluent®,  will  be  used  for  these  simulations  and  Gridgen® 
will  be  used  to  create  the  grids. 

The  research  presented  requires  low  velocities  at  ~2.0  m/s  (near- stationary  flow) 
and  v  =  =  {15.2,  30.4)  m/s  in  order  to  scope  the  trade  space  and  compare  against 

experimental  data  displayed  in  Figure  12,  Figure  13,  Figure  14,  and  Figure  15.  These 
velocities  correspond  to  incompressible  flows.  A  validated  incompressible  flow  solver 
would  take  a  significant  amount  of  time  to  create,  much  more  than  is  reasonable  for  the 
purposes  of  this  research  effort.  This  is  the  main  reason  for  employing  Fluent®.  Fluent® 
is  a  commercial  software  package  that  can  solve  2-D  and  3-D  fluid  flow  simulations.  It 
can  handle  a  wide  variety  of  flow  conditions,  such  as  compressible  and  incompressible 
flows.  An  implicit  incompressible  method  of  an  unsteady  time- accurate  solution  will  be 
used.  Within  each  time  step,  sub -iterations  may  be  performed  to  reduce  the  residual.  A 
maximum  of  20  sub-iterations  or  a  tolerance  of  10‘6  for  the  residual  will  be  used  before 
the  solver  continues  to  the  next  time  step.  Simulation  of  the  DBD  operation  will  be 
performed  with  a  set  of  User  Defined  Functions  (UDF’s)  implemented  in  Fluent®. 


NACA  0009 


Figure  12 

Lift  Coefficient  vs  Angle  of  Attack  for 
Rec=180k  without  and  with  Actuator 
Operating  [7] 


NACA  0009 


Figure  13 
Drag  Polar  for 

Rec=180k  without  and  with  Actuator 
Operating  [7] 


NACA  0009 


Figure  14 

Lift  Coefficient  vs  Angle  of  Attack  for 
Rec=360k  without  and  with  Actuator 
Operating  [7] 


NACA  0009 


Figure  15 
Drag  Polar  for 

Rec=360k  without  and  with  Actuator 
Operating  [7] 


First,  before  going  into  the  simulation  of  a  DBD  using  UDF’s,  it  will  be  necessary 
to  establish  validation  of  the  test  cases.  The  first  set  of  test  cases  will  be  a  simple  flat 


1-11 


plate.  The  second  set  of  test  cases  will  be  flow  over  a  NACA  0009  airfoil  with  the 
specifications  from  Figure  16. 


40  60 

x/c  (%) 


too 


Airfoil  Type  NACA  0009 
Chord  20.2  cm 

Span  41.7  cm 

Velocity  15.2,  30.4  m/s 

Rec  0.18x10®,  0.36x10® 

p *  0.993  kg/m3 

i/*  2.025  x  10~3  m2/s 


rBased  on  7000  ft  altitude. 


Figure  16 

NACA0009  Airfoil  Test  Parameters  [7] 


The  flat  plate  will  be  validated  by  subjecting  the  grid  to  both  of  the  above 
Reynolds’  cord  numbers  for  the  airfoil,  where  the  Reynolds’  number  will  be  assumed  to 
be  the  same  along  a  flat  plate  as  along  an  airfoil,  Ref  =  Rev  =  |0.18xl06,  0.36xl06  j .  The 

boundary  layer  profile  will  be  compared  to  the  analytic  Blasius  differential  equations 
solution  for  validation 

The  airfoil  will  be  validated  by  subjecting  the  grid  to  a  set  of  simulations  with 
v  =  E/0O  =  {15.2,  30.4]  m/s,  corresponding  to  Rec  =|0.18xl06,  0.36xl06] ,  at  angles  of 

attack  spanning  -16°  to  +16°  in  1°  increments.  The  coefficient  of  lift,  Cl,  and  the 
coefficient  of  drag,  Cd,  will  then  be  compared  to  experimentally  known  data  obtained 
from  Selig  [15]  for  validation. 

There  are  many  different  types  of  UDF’s  for  Fluent®.  Source  Term  UDF’s  will 
be  used  to  calculate  the  simulated  addition  of  momentum  and  thermal  energy  to  the  flow. 
A  Define  on  Demand  UDF  will  be  used  to  spatially  distribute  the  momentum  and  thermal 
energy  addition  to  the  flow.  Hie  UDF’s  are  written  in  the  C  programming  language. 
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These  subroutines  will  allow  a  macroscopic  simulation  of  the  DBD  and  are  the  purpose 
of  this  research  effort.  More  detail  on  these  subroutines  and  their  design  can  be  found  in 
the  Simulation  Setup  Section  and  Appendix  B. 

Corke  et  al.  [7]  present  the  main  mechanism  contributing  to  the  boundary  layer 
flow  as  force  created  via  ion- neutral  collisions  from  the  positive  ions  accelerating  in  the 
electric  field.  A  simplified  equation  for  the  pressure  term  coupling  to  the  neutral  gas  flow 
is  given  by  Corke  et  aL  [7]  in  Equation  (1). 


BE=|e0V£2 


(1) 


The  unknown  variable  of  this  equation  is  the  electric  field,  which  includes  not 
only  the  induced  field  from  the  potential  between  the  electrodes,  but  also  includes  the 
field  from  the  plasma  as  well.  As  a  result,  we  refer  to  an  updated  version  of  Equation  (1) 
in  Equation  (2)  from  Corke  et  al.  [10],  which  includes  the  charge  density  of  the  plasma 
estimated  from  an  electrostatic  view,  but  is  still  an  intuitive  approximation 


/;=pc£=- 


vV; 


cp  E 


(2) 


In  contrast  with  the  previous  two  equations,  Boeuf  and  Pitchford  [16]  use  an 
approach  relating  “...the  force  per  unit  volume  acting  on  the  gas  molecules...”  with  the 
number  of  ions,  p,  the  electric  field,  E ,  the  ion  current  density,  ji,  and  the  ion  mobility, 
[Li,  in  Equation  (3). 


(3) 
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This  equation  assumes  that  the  force  is  primarily  transferred  in  a  non-neutral  region 
during  an  ion  and  neutral  particle  collisions  where  the  ion  number  density  is  much  greater 
than  the  electron  number  density. 

The  calculations  for  solving  the  electric  field  variable,  E ,  in  Equations  (1),  (2), 
and  (3)  are  quite  involved  and  require  that  the  time  steps  taken  be  small  with  respect  to 
the  time  of  one  wave  cycle.  The  frequencies  used  generally  reside  between  1  kHz  and  10 
kHz.  Frequencies  such  as  these  require  a  significant  amount  of  computational  time  to 
arrive  at  a  valid  simulation.  Therefore,  a  macroscopic  approach  of  the  average  effects 
due  to  the  DBD  over  several  wave  cycles  is  desirable. 

The  purpose  of  this  research  is  to  examine  a  macroscopic  view  (force  and  power 
addition)  of  a  DBD  in  operation.  The  Source  Term  UDF  will  add  momentum  and  energy 
(time  derivatives  of  force  and  power  respectively)  to  the  flow  in  an  attempt  to  model  the 
behavior  of  a  DBD  without  solving  a  complex  and  calculation  intensive  equation  A 
“weighting”  function  will  assign  values  to  each  cell  in  order  to  distribute  the  momentum 
and  energy  addition  over  a  particular  spatial  extent.  Further,  the  time  for  one  period  of 
the  AC  waveform  at  a  kHz  frequency  is  still  much  smaller  than  the  anticipated  time  step 
associated  with  the  flow  simulation.  As  such,  a  temporal  average  of  the  force  over 
several  cycles  will  be  used  and  implemented  via  the  weighting  function  There  are 
currently  two  different  views  as  to  how  to  spatially  distribute  the  source  terms  imparted 
from  the  DBD  on  the  simulation  grid  as  well  as  their  values. 
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Boeuf  and  Pitchford  Impulse  Density 

The  first  distribution  is  to  employ  the  source  term  in  an  extremely  localized  set  of 
cells  according  to  Boeuf  and  Pitchford’s  impulse  simulations  [16].  Their  simulations 
were  performed  with  Nitrogen  gas  at  Standard  Temperature  and  Pressure  (STP)  and 
account  for  secondary  ionization.  The  model  used  does  not  employ  a  neutral  gas  flow 
solver;  instead,  it  is  an  ionized  gas  solver  code  adapted  from  their  extensive  experience 
with  plasma  display  panels. 

The  device  geometry  illustrated  in  Figure  17  has  length  scales  that  are  small,  200 
pm  by  800  pm,  with  equally  small  cell  sizes  (not  illustrated  in  Figure  17)  of  2  pm  on  a 
side.  They  state  “. .  .the  average  force  per  unit  volume. .  .will  be  in  the  10'  -  10  N  m" 
range”.  The  median  of  this  range  was  estimated  to  be  consistent  with  the  results  of  Font 

[8],  103  N/m3 . 
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Figure  17 

Boeuf  and  Pitchford  Simulation  Geometry 

Boeuf  and  Pitchford  also  report  the  contours  of  the  impulse  density,  F  ■  At , 

around  a  DBD  as  shown  in  Figure  18  and  Figure  19.  Their  simulations  used  a  single 

square  wave.  If  they  used  a  1  kHz  driving  square  wave,  the  impulse  density  profile 

would  be  multiplied  by  1000  (=1  kHz)  to  give  a  force  density. 
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Position  (urn) 

Figure  18 

Boeuf  and  Pitchford  [16] 
X-Component  Impulse  Density  Weighting 
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Figure  19 

Boeuf  and  Pitchford  [16] 

Y- Component  Impulse  Density  Weighting 


Finally,  Boeuf  and  Pitchford  estimate  that  the  maximum  increment  to  the  velocity 
magnitude  of  the  fluid  is  directly  related  to  the  X- momentum  impulse  weight  by 
av  =  -i-J*/  dt ;  where  the  fluid  density  is  p  =  1.2  kg/m3 ,  /  is  the  impulse  density,  and  dt 

is  time.  Using  Figure  18  will  produce  a  wall  jet  with  a  velocity  between  5  m/s  and  10  m/s 
at  a  height  of  1  Op  m  off  the  surface  of  the  flat  plate.  The  boundary  layer  profile 
associated  with  this  effect  is  depicted  in  Figure  20  and  Figure  21. 


1-16 


Y  (m)  Y  (m) 


0.00600 


Figure  20 

Boeuf  and  Pitchford  Estimation  of  Wall-Jet  Velocity 


Figure  21 

Boeuf  and  Pitchford  Estimation  of  Wall-Jet  Velocity  (Close-up) 
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Roy  and  Gaitonde  Force  Density 

The  second  method  of  distributing  the  force  and  thermal  addition  densities  is  by 
Roy  and  Gaitonde  [17],  who  show  that  the  force  density  at  the  DBD  is  on  the  order  of 
103  pN/cm3  ( 1 03  N/m3 ) .  The  cm3  unit  volume  was  confirmed  with  the  authors  even 

though  it  was  not  stated  specifically  in  their  publication.  Further  discussions  with  the 
authors  revealed  that  the  input  power  to  the  system  per  unit  length  was  approximately 
7  W/m  with  negligible  amounts  of  this  power  going  towards  thermal  heating.  Boundary 
layer  velocity  profiles  reaching  2.5  m/s  are  shown  in  Figure  25  and  are  further  discussed 
in  the  Results  and  Conclusions  Sections.  Their  simulations  were  performed  using 
Helium  at  STP  without  secondary  emission.  The  use  of  Helium  instead  of  a  diatomic 
molecule  may  have  a  significant  impact  on  the  model’s  performance  and  may  not  mimic 
atmospheric  gas  effects  correctly  as  a  result. 

The  Roy  and  Gaitonde  distribution  of  force  density  is  several  orders  of  magnitude 
larger  than  the  Boeuf  and  Pitchford  distribution  The  dimensions  for  the  computational 
volume  under  the  DBD  influence  are  0.5  cm  high  by  3  cm  wide  by  1  m  deep  for  the  Roy 
and  Gaitonde  model,  compared  to  200  pm  (0.02  cm)  high  by  800  pm  (0.08  cm)  wide  by 
1  m  deep  for  the  Boeuf  and  Pitchford  model.  Further,  the  widths  of  the  electrodes  for 
each  case  vary  in  the  same  respect.  The  Roy  and  Gaitonde  geometry  uses  electrodes  that 
are  1.2  cm  wide,  while  the  Boeuf  and  Pitchford  geometry  uses  electrodes  that  are  100  pm 
(0.01  cm)  wide  for  the  exposed  electrode  and  an  800  pm  (0.08  cm)  wide  buried  electrode 
that  spans  the  entire  simulation  space.  The  disparity  in  size  of  the  simulation  space  as 
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well  as  the  electrodes  is  assumed  to  have  an  effect  on  the  performance  of  the  two 
systems. 

Data  obtained  from  Roy  and  Gaitonde  show  their  simulation  results  for  the 
induced  force  densities  in  Figure  22,  Figure  23,  and  Figure  24.  The  bold  red  lines  in 
these  figures  represent  the  electrodes  of  the  DBD.  As  was  already  reported,  these 
electrodes  are  1.2cm  wide,  two  times  wider  than  the  electrodes  that  were  used  in 
experiments  run  by  Post  and  Corke  [4]  using  a  NACA  0009  airfoil.  Length  scales  are 
compared  and  contrasted  later  in  the  validation  section  Roy  has  hypothesized  that  the 
long  length  of  the  lower  electrode  allows  a  charge  buildup  on  the  dielectric  surface  that 
creates  the  negative  force  depicted  in  Figure  22  and  Figure  23.  Figure  25  depicts  the  wall 
jet  velocity  profile  at  2mm  increments,  starting  from  2mm  upstream  of  the  DBD 
electrode  juncture. 


Figure  22 

Macroscopic  View  of  X- momentum  Force  Field  (N/m3)  [17] 
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Figure  23 

O 

X-momentum  Force  Field  (N/m  )  [17] 


cm 


Figure  24 

Y-momentum  Force  Field  (N/m  )  [17] 
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Computed  Streamwise  Velocity  Induced  in  a  Quiescent  Helium  Gas  [17] 


From  the  three  sources  examined:  Font  [8],  Boeuf  and  Pitchford  [16],  and  Roy 
and  Gaitonde  [17] ;  the  force  densities  appear  to  be  set  around  103  N/m3.  Therefore,  it  is 
suggested  that  the  momentum  source  term  defined  in  units  of  force  per  unit  volume 
should  be  near  103  N/m3  to  produce  a  wall  jet  of  ~5  m/s  in  a  near- stationary  flow. 

Data  Set  Test  Plan 

The  first  set  of  test  cases  will  simulate  flow  response  using  the  provided  force 
densities  from  Boeuf  and  Pitchford  [16]  and  Roy  and  Gaitonde  [17]  on  the  flat  plate 
grids.  As  each  scheme  requires  different  length  scales  for  the  grid  cell  sizes,  two 
different  yet  fundamentally  similar  grids  will  be  required.  The  geometry  of  the  grids  for 
each  of  these  schemes  is  explained  in  detail  in  the  Validation  section. 

After  the  initial  force  density  profiles  are  both  complete,  several  more  simulations 
will  be  run  in  order  to  explore  the  force  density  magnitude  in  each  scheme  to  sufficiently 
induce  a  5  m/s  wall  jet  in  a  flow  of  2  m/s  .  A  simple  analysis  of  heating  and  momentum 
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addition  is  performed  in  the  Simulation  Setup  section.  Until  this  point,  little  to  no 
thermal  energy  will  have  been  put  into  the  flow. 

Next,  thermal  energy  input  will  be  increased  into  the  flow,  showing  the  effects  of 
thermal  energy  on  the  wall  jet  characteristics.  A  set  of  simulations  will  also  be  run  where 
the  percentage  of  the  components  of  thermal  and  momentum  addition  are  kept  constant 
while  the  total  power  is  increased. 

Finally,  a  NACA  0009  airfoil  will  be  simulated  with  and  without  the  additionof 
the  DBD  simulation  source  terms.  Simulated  lift  and  drag  characteristics  will  be 
compared  to  experimental  data  in  Figure  12  thru  Figure  15  as  reported  by  Corke  [7]. 

Expectations 

This  research  should  establish  a  macroscopic  (force  and  power  addition) 
computational  simulation  of  a  DBD’s  momentum  and  thermal  energy  transfer  to  a  flow. 
Boundary  layer  velocity  profiles  will  be  obtained  and  compared  with  the  experimental 
data  reported  by  Corke  [7]  and  Newcamp  [18].  The  question  of  “How  much  momentum 
and  thermal  energy  are  imparted  to  the  flow  and  in  what  fashion?”  is  to  be  answered. 


1-22 


II.  Simulation  Setup 


Several  pieces  of  code  needed  to  be  developed  in  order  to  allow  the  commercial 
software,  Fluent®,  to  accurately  simulate  a  DBD.  First,  there  is  a  subroutine  that 
Fluent®  calls  once,  just  prior  to  starting  a  simulation.  This  subroutine  stores  a  “weight” 
value  for  each  cell  in  the  nearby  vicinity  of  the  DBD  location  so  that  the  three  main 
subroutines  can  quickly  have  access  to  a  set  of  normalized  “weighting”  data.  The  three 
main  subroutines  are  called  every  time  an  iteration  computes  cell  data.  The  subroutines 
involve  the  local  addition  of  thermal  energy,  x- momentum,  and  y- momentum  into  the 
system. 

Each  time  step  was  0.001  seconds.  This  is  approximately  the  amount  of  time 
information  in  the  fluid  takes  to  travel  the  length  of  the  flat  plate  or  airfoil,  0.202  meters, 
computed  by  dividing  the  cord  length  by  the  speed  of  sound.  Smaller  time  steps  may 
increase  accuracy  of  the  simulation,  but  they  always  increase  the  computational  time.  A 
total  of  1000  time  steps  were  completed  for  a  total  time  of  1  second.  Up  to  20  sub- 
iterations  are  performed  for  every  time  step,  with  convergence  of  the  residual  being 
monitored.  If  the  residual  became  less  than  10"6  or  all  20  sub -iterations  were  completed, 
then  the  simulation  proceeded  to  the  next  time  step.  Monitoring  the  residual  has  the 
practicality  for  aiding  in  a  more  accurate  solution  using  an  iterative  approach,  while  still 
time  stepping  so  as  to  examine  unsteady  phenomenon. 

Cell  Weighting  Subroutine 

The  subroutine  that  weights  cells  within  the  DBD  vicinity  is  extremely  important. 
The  spatial  extent  of  the  DBD  changes  depending  upon  which  force  density  profile  is 


chosen;  Boeuf  and  Pitchford  or  Roy  and  Gaitonde.  Flexible  code  is  required  to  adapt  to 
each  force  density  profile  and  to  span  over  several  grid  cells,  no  matter  their  size  and 
shape.  The  DEFINE_ON_DEMAND(cell_weight_on_demand)  subroutine  exercises  this 
function 

The  UDF  C  code  is  written  in  several  subroutines,  all  of  which  are  detailed  in 
Appendix  B.  First,  the  given  location  of  the  DBD  is  taken  from  the  code  and  the  tangent 
vector  to  the  surface  at  that  location  is  found.  A  coordinate  transformation  is  then 
employed  so  that  this  vector  forms  the  new  x-axis.  Each  transformed  point  is  tested  to 
see  if  it  lies  within  an  estimated  influence  boundary  of  the  DBD,  as  predetermined  by  the 
force  density  profile  that  was  chosen  for  modeling.  The  power  density  for  thermal 
addition  uses  the  force  density  profile  to  distribute  energy,  as  there  is  no  available  data 
for  how  the  thermal  addition  is  distributed.  If  the  point  does  lie  within  an  estimated 
influence  boundary,  then  a  series  of  calculations  occur  to  give  the  cell  a  weight.  Many 
issues  are  taken  into  account,  such  as  cells  straddling  one  or  more  boundaries  with 
different  equations  describing  each  boundary  area.  This  was  one  of  the  most  complex 
parts  of  this  thesis  to  design.  Without  it,  the  developed  code  would  be  too  rigid  for  any 
follow-on  work. 

The  cell  weighting  equation  set  can  produce  force  density  profiles  such  as  those 
displayed  in  Figure  26  and  Figure  27.  Distances  are  in  meters,  with  the  DBD  location  set 
at  75%  cord  length,  or  x  =  0.1515  meters  for  a  cord  length  of  0.202  meters.  The  cell 
weighting  equation  set  can  be  found  in  Appendix  B  in  the  Subroutine:  weight_funct 
section  The  equations  that  represent  the  local  weight  were  created  using  the  force 
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density  plots  from  Roy  and  Gaitonde  and  impulse  density  plots  from  Boeuf  and 


Pitchford.  Equation  (4)  was  used  to  create  a  Boeuf  and  Pitchford  force  density  profile. 
Similarly,  Equation  (5)  was  used  to  create  a  Roy  and  Gaitonde  force  density  profile.  In 
each  of  these  equations,  the  location  of  the  DBD’s  center  is  defined  as  (x,y)  =  (0,0).  All 
cells  outside  of  these  ranges  are  given  a  weight  of  zero. 

lOOOOexp [-50000 |x| - 80000 |y|]|  J-0.0001m<x  <  +0.0000m 

lOOOOexp  [-100  |x|-  80000  |y|]  °'  {+0.0000m  <  x  < +0.0006m  (4) 

-0.001m  <  y  <  +0. 00005m 


163exp[-5  x  -5|y|] 


3200*  1 0-275'+275W 

100*1  o“240ltl“500l-vl 


y  for  ( 


-0.0135m  <  x  <  -0.0045m 
-0.0045m  <  x  <  +0. 0040m 
+0. 0040m  <  x  <  +0.01 65m 


-0.001m  <  y  <  +0.005m 


(5) 


The  values  in  Equations  (4)  and  (5)  were  found  using  a  parametric  approach  to 
model  the  behavior  of  a  DBD.  The  ranges  were  derived  from  Figure  18  and  Figure  19  for 
Equation  (4),  and  Figure  23  and  Figure  24  for  Equation  (5).  As  raw  data  was  not 
obtained  from  either  Boeuf  and  Pitchford  or  Roy  and  Gaitonde,  a  direct  comparison 
between  the  original  force  density  profiles  and  the  fitted  force  density  profile  is  not 
available. 

The  cell  weighting  code  also  takes  into  account  the  shape  of  the  airfoil  at  the 
location  of  the  DBD.  This  added  computational  effect  can  best  be  seen  when  comparing 
the  slopes  of  the  gradient  contour  boundaries  between  Figure  26  and  Figure  27. 
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Figure  26 

Example  of  Flat  Plate  Cell  Force  Density  Profile 
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Figure  27 

Example  of  NACA0009  Airfoil  Cell  Force  Density  Profile 


Once  the  weight  of  all  of  the  cells  has  been  determined,  the  total  weight  is 
computed  and  used  to  normalize  all  of  the  cells’  weights.  Cells  outside  the  weighted 
boundary  are  given  a  weight  of  zero.  Finally,  the  cell  weight  is  stored  at  the 
corresponding  cell’s  center  where  Fluent®  can  call  up  the  data  when  running  the  other 
three  main  subroutines. 
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Source  Simulation  of  a  DBD 


To  simulate  a  DBD,  three  subroutines  were  created  that  add  thermal  energy 
(temp_source),  x- momentum  (x_momentum_source),  and  y- momentum 
(y_momentum_source)  to  the  flow.  None  of  these  subroutines  add  mass  to  the  flow,  an 
approach  which  is  typical  for  a  wall  jet  simulation,  but  not  for  a  DBD  simulation. 

The  DBD  model  has  to  account  for  certain  parameters  that  are  controlled  during 
physical  experimentation.  The  input  power  in  Watts  can  not  be  exceeded  in  the 
simulation,  which  is  specified  in  Watts  per  unit  span  length  of  1  meter  for  the  2D 
simulations.  Because  all  three  of  the  source  additions  must  sum  to  the  input  power,  they 
must  be  divided  by  the  following  method.  A  percentage  of  the  input  power  is  given  to 
thermal  energy,  and  the  remaining  amount  necessary  to  give  100%  total  power  usage  was 
given  to  momentum,  as  seen  in  Equation  (6).  The  power  per  length  in  W/m  for  each  of 
these  portions  is  defined  as  Thermal  Power  and  Momentum  Power,  respectively.  Each 
momentum  component  is  then  further  divided  from  the  power  given  to  momentum, 
where  the  total  percentage  of  the  x  and  y  components  of  momentum  addition  are  to  equal 
100%,  as  seen  in  Equation  (6).  These  percentages,  multiplied  by  the  cell  weight 
calculated  in  the  cell  weighting  subroutine,  give  the  total  amount  of  power  delivered  to  a 
particular  cell  in  terms  of  power  density  in  units  of  W/m  and  directional  force  density  in 
units  of  N/m  .  Watts  and  Newtons  are  the  units  of  the  time  derivative  for  energy  and 
momentum  respectively.  A  derivation  for  the  transformation  of  power  into  power  density 
and  force  density  for  thermal  energy  addition  and  momentum  addition,  respectively, 


follows. 


Thermal  Energy  %  +  Momentum  %  =  100% 
X-Momentum  %  +  Y-Momentum  %  =  100% 


(6) 


Thermal  Energy  Power  +  Momentum  Power  =  Total  Input  Power 
X-Momentum  Power  +  Y-Momentum  Power  =  Momentum  Power 

If  the  thermal  energy  addition  has  a  considerable  effect,  then  a  radial  expansion  of 
the  flow  will  be  expected  where  the  free  stream  flow  velocity  is  stationary,  as  in  Figure  3. 
This  is  not  what  is  depicted  by  Digital  Particle  Image  Velocimeter  (DPIV)  measurements 
from  Figure  3.  Instead,  the  flow  is  accelerated  to  the  right,  which  is  consistent  with 
momentum  addition  in  the  positive  x- direction.  There  is  also  a  component  of  momentum 
addition  in  the  y-direction  as  the  flow  is  drawn  towards  the  DBD. 

Extreme  care  must  be  taken  when  implementing  and  reporting  source  terms  in 
Fluent®.  The  documentation  and  examples  supplied  by  the  company  and  found  online 
are  misleading.  The  units  for  the  thermal  energy  UDF  source  and  the  momentum  UDF 
source  may  be  thought  to  be  J/m3  and  kg  •  (m/s)/m3 ,  respectively.  In  actuality,  the 
correct  units  are  time  derivatives  of  these  units  as  W/m  and  N/m  ,  respectively. 

The  correct  method  for  thermal  energy  addition  to  the  flow  is  simple.  Take  the 
power  in  Watts  to  be  delivered  as  thermal  energy  addition  and  divide  by  the  cell  volume 
to  produce  the  required  input  for  thermal  energy  addition  with  units  of  W/m .  Equation 
(8)  shows  this  in  equation  form.  Due  to  a  lack  of  experimental  and  simulated  power 
density  distribution  data,  it  was  assumed  that  the  power  density  distribution  was  similar 
to  the  force  density  distribution.  This  may  or  may  not  be  accurate,  and  an  examination  of 


the  effect  of  power  density  to  the  system  is  performed  in  the  Results  and  Conclusions 
sections. 

_  __  .  Cell  Weight  *Thermal  Energy  %*Total  Power 

Power  Density= - - - - - - - - - -  (8) 

Cell.  Volume 

A  procedure  for  momentum  addition  was  more  difficult  to  establish.  A  one 
dimensional  case  will  be  examined  for  equation  development.  A  two  dimensional  case 
would  require  adding  a  similar  set  of  equations  for  the  added  dimension.  Also,  because 
there  is  a  source  term  subroutine  for  each  momentum  addition  component,  only  one 
velocity  component  needs  to  change  per  momentum  subroutine.  The  kinetic  energy  of 
the  flow,  KEo,  was  first  solved  for  using 

KE0  =  jtnu02  (9) 

where  uo  is  a  one  dimension  velocity  component,  m  is  the  mass  within  the  cell  computed 
from  m  -  p  Cell_Volume,  and  p  is  the  fluid  density  in  the  cell.  Once  the  fluid’s  kinetic 

energy  is  found,  the  equation  used  to  derive  incremental  energy  from  power  is  computed 
in  Equation  (11).  To  allow  for  the  flow  direction  to  have  effect,  the  change  in  kinetic 
energy,  aKE  ,  is  defined  as  positive  for  forces  directed  from  left  to  right,  and  negative  for 
forces  directed  from  right  to  left.  This  is  due  to  the  cell  weight  introducing  a  positive  or 
negative  sign  as  a  force  density  profile  requires.  The  total  power  is  a  user  input 
parameter  for  the  simulations  and  must  be  positive.  In  order  to  check  that  the  program  is 
performing  properly,  the  calculation  in  Equation  (10)  may  be  performed. 

Total_Power  =  £  (|AKEn  |  At) 
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(10) 


aKE  =  Power*At 


where  Power  =  Cell_Weight*Momentum_%*Total_Power 


(11) 


The  aKE  is  then  added  to  the  old  kinetic  energy,  KEo,  to  give  a  new  kinetic  energy,  KE] . 

(12) 


KEj  =jm«02  +aKE 


Equation  (13)  is  used  to  preserve  direction  of  a  velocity  component.  The  inverse  of 
Equation (13)  is  Equation  (14),  where  a  =u-  \u\ . 


u  =>u-\u\ 


(13) 


lfa<0  M=>-J|a| 
If  a  >  0  w  =>  +-Ja 


(14) 


Briefly,  two  cases  will  prove  Equation  (13)  and  (14).  If  u  =  1 ,  then  u  =1-1  =  1  and  the 
sign  is  preserved.  If  «  =  (-l),then  u2  =  (-1)- 1(-1)|  =  (  -1)  -|l|  =  -l  and  the  sign  is  again 

preserved.  Further,  if  a  =1 ,  then  u  =  -K/l  =1  and  the  sign  is  preserved.  Similarly,  if 
a  =  (-1) ,  then  u  =  (-1)  =  -~J\  -  (-1)  and  the  sign  is  again  preserved. 


The  new  kinetic  energy  of  the  flow  can  again  be  described  using  mass  and  velocity. 

KE,  =-(■ mu 2 

Therefore,  Equation  (13)  is  substituted  into  Equation  (15)  to  produce 

KEt  =\mu2  =-jm.[ui  ■  | ) 


(15) 


(16) 


Solving  for  the  velocity  component  yields 


2  -  KE, 


u\  \u\  =" 


(17) 


m 
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Substituting  in  Equation  (12)  for  KEj 


u\ '  Mi  =' 


2'(lm(uo  •|mo|)+aKE) 


Simplifying  produces 


ux  ■  Mj  |  ={u0  ■  |^o  |) 


2  aKE 


If  the  quantity  in  Equation  (19)  is  negative,  then  an  additional  absolute  value  and 
negative  sign  must  be  implemented  to  keep  the  velocity  sign  consistent  as  was  shown 
with  Equation  (14).  This  is  shown  with  Equation  (20).  If  the  new  flow  velocity  is 
positive,  or  flows  to  the  right,  then  Equation  (21)  would  be  used. 


ui  =  ~.  (mo'K|)- 


2  ■  (Cell_Weight*Momentum_%*Total_Power ) 
p  CellJVolume 


u,  =  +. 


( uo ■  bl)- 


2  ■  (Cell_Weight*Momentum_%*Total_Power ) 
p  Cell_Volume 


The  difference  between  the  new  and  the  old  velocity  is  then  divided  by  the  time  step  to 
yield  acceleration. 


Wj  u0 


The  acceleration  multiplied  by  density  gives  the  source  addition  term  in  units  of  N/m3. 
Fluent®  requires  that  a  source  be  in  terms  of  a  change  per  volume,  which  yields 


m- a  u,-i i 

- =  p -a  =p  ■— — 

vol  At 
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Equation  (24)  shows  this  entire  procedure  in  one  form  with  the  assumption  that  the 
kinetic  energy  increase,  aKE  ,  and  the  initial  starting  velocity  are  >  0 . 

j  2  •  (Cell_Weight*Momentum_%*Total_Power)  •  At  2 

\  0  p  Cell_Volume  0  (24) 

Force  Density=p  ^ - — - 2 - 

As  a  check  for  the  Roy  and  Gaitonde  force  density  profile,  the  velocity  increase 
was  calculated  using  a  time  step  of  0.001  seconds  and  2%  momentum  from  5  W/m  ,  or 
0.10  W/m  ;  a  setting  that  yields  a  force  density  of  1800  N/m  using  Equation  (24)  with  a 
stationary  flow,  uQ  =  0  m/s.  This  is  consistent  with  the  Roy  and  Gaitonde  force  density 
profile.  To  compute  this  increase,  the  volume  was  computed  from  the  dimensions  of  the 
weighted  area  of  the  model.  The  dimensions  are  0.5cm  x  l.Ocmx  1 00cm (l  meter) ,  to 
give  a  volume  of  50xl0“6m3  =  50cm3 .  The  density  at  STP  is  p  =1.225  kg/m3 ,  and  the 
mass  in  the  volume  is  found  by  m  =  p  •  vol .  Assuming  a  stationary  flow,  uQ=  0  m/s, 
using  Equation  (19)  and  the  positive  part  of  Equation  (14),  the  result  yields  an  estimated 
increase  of  1.8  m/s. 


As  a  check  for  the  Boeuf  and  Pitchford  force  density  profile,  the  velocity  increase 
was  calculated  using  a  time  step  of  0.001  seconds  and  0.00027  W/m;  a  setting  that  will 
later  prove  in  the  Results  section  to  yield  a  force  density  consistent  with  the  Boeuf  and 
Pitchford  force  density  profile.  The  dimensions  for  the  Boeuf  and  Pitchford  geometry  are 
150|lmx  800p,mx  1 00cm (l  meter) ,  to  give  a  volume  of  120x10  9m3  =  120mm3 .  The 
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same  density,  stationary  flow  equation,  and  Equations  (24)  and  (25)  were  used  to  yield  a 
force  density  of  1900  N/m  and  giving  a  result  of  1.9  m/s  for  the  wall  jet  velocity, 
respectively. 


Further,  when  calculating  using  a  moving  fluid  the  equation  requires  the 
additional  variable  of  initial  flow  velocity,  u0 .  Using  the  same  assumptions  of  positive 
flow  velocity  and  positive  aKE  the  equation  now  becomes 


u  —  .  u0~  +  ■ 


2-aKE 


2  ■  Pwr  ■  At 


m 


p  vol 


(26) 


If  the  increment  of  energy  to  the  cell  is  fixed  by  a  force  density  profile,  then  a 
relationship  between  the  new  velocity  and  the  initial  may  be  graphed,  as  is  done  in  Figure 
28  and  Figure  29. 


— 1 ■—  Roy  &  Gaitonde  — —  Boeuf  &  Pitchford 


Figure  28 

Estimated  Wall  Jet  Peak  Velocity  Magnitude  (m/s) 
Compared  to  the  Free  Stream  Velocity  (m/s) 


n-n 


—> ■—  Roy  &  Gaitonde  Difference  — *—  Boeuf  &  Pitchford  Difference 


Figure  29 

Difference  between  the  Estimated  Wall  Jet  Peak  Velocity  Magnitude 
and  the  Free  Stream  Velocity  (m/s) 

Compared  to  the  Free  Stream  Velocity  (m/s) 


As  with  the  experiments,  as  the  free  stream  velocity  is  increased,  the  wall  jet 
becomes  less  notable.  This  is  because  the  increment  in  energy  delivered  to  the  flow 
remains  the  same,  but  the  fluid’s  kinetic  energy  increases  as  the  free  stream  velocity 
increases.  Eventually,  the  energy  increment  becomes  insignificant  when  compared  to  the 
kinetic  energy  of  the  fast  moving  free  stream  velocity  and  little  to  no  change  in  the 
velocity  profile  is  seen. 

Now  that  the  momentum  and  energy  addition  subroutines  have  been  explained 
and  their  effects  on  the  flow  estimated,  it  is  necessary  to  describe  the  term  ds[eqn].  In  all 
three  of  the  subroutines,  the  term  ds[eqn]  is  set  to  0.0.  This  term  is  the  derivative  of  the 
source  term  if  known  or  0.0  if  unknown.  Because  the  derivatives  of  the  source  terms  are 
complex  and  dependent  upon  the  model,  Fluent®  was  employed  to  explicitly  solve  for 
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the  derivative  in  its  attempt  to  derive  a  better  and  more  stable  solution  for  its  implicit 
solver.  A  sample  momentum  source  term  that  the  Fluent®  manual  supplied  was  run  with 
ds[eqn]  set  equal  to  the  derivative  of  a  source  and  ds[eqn]=0.0  in  a  second  case.  Both 
cases  converged  in  the  same  number  of  iterations  and  had  the  exact  same  results.  The 
clock  time  difference  was  measured  as  insignificant  for  this  test  run  set  and  is  not 
expected  to  have  any  significant  impact  on  simulation  run  times. 

Finally,  to  ensure  that  the  additional  subroutines  written  were  all  functioning 
correctly,  a  set  of  verification  runs  outside  of  Fluent®  were  performed.  Further  detail  of 
each  subroutine  and  its  workings  are  discussed  in  Appendix  B.  The  verification  for  each 
subroutine  is  also  discussed  there  as  well. 
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III.  Validation 


Calculation  for  Boundary  Layer  Thickness  and  Flat  Plate  Grid  Geometry 

For  flow  along  a  flat  plate  as  seen  in  Figure  30,  the  boundary  conditions  are  no¬ 
slip  at  the  wall,  u{ 0)  =  0 ,  and  a  smooth  transition  of  the  flow  to  the  free  stream  velocity, 

u  (699 )  =  6%  and  ^  =  0  .  If  the  flow  is  laminar,  these  boundary  layer  velocities  can 

ay  . 

V=099 

be  approximated  with  a  second  order  polynomial  approximation  in  Equation  (27),  given 
by  White  [1 1 :222],  equation  4-11.  The  Blasius  formula  for  determining  where  the 
boundary  layer  thickness  is  99%  of  the  free  stream  flow  is  given  in  Equation  (28),  from 
White  [1 1 :223],  equation  4- 14;  where  p  is  density  with  units  of  kg/m  ,  v  =  U„  is  the  free 

stream  velocity  with  units  of  m/s,  ji  is  viscosity  with  units  of  Ns/mr,  and  x  is  the  distance 
in  meters  from  the  leading  edge  of  the  flat  plate.  Equation  (29)  is  the  Reynold’s  number 
at  a  point  x  along  the  flat  plate  from  the  leading  edge. 


u{y)~uo „ 


(27) 


(28) 
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Boundary  Layer  Velocity  Profile  for  a  Flat  Plate 
The  goal  for  validation  is  to  model  the  DBD  operation  at  near  stationary  flow 
velocities  at  STP  and  with  the  same  setup  as  the  NACA0009  airfoil  experiments.  From 
the  table  in  Figure  16,  we  can  extract  p  =  0.993  kg/m3 ,  v  =  =  {15.2,  30.4]  m/s ,  the 

cord  length  c  =  0.202m,  and  Rec  =|0.18xl06,  0.36xl06]  for  the  NACA0009  airfoil 

experiments  reported  by  Corke  [7].  Because  the  viscosity  at  7000  feet  altitude  was  not 
given,  the  Rex  equation  was  inverted  for  p,  giving  p  =  1 .69x10  5  N-  s/m2 .  This 


corresponds  to  kinematic  viscosity,  u  =  —  =  1.70x10  5  m2/s  .  The  temperature  was  not 

P 

published  and  will  be  assumed  to  be  288.15  °K,  with  a  pressure  of  78669  Pa  based  on  a 


7000  foot  (2133.5  meters)  altitude,  derived  from  the  equation  p  =  pQ  exp 


~g 

R*T Z 


where  p0  =101325  Pa,  R  =  287  J/kg ° K ,  g  =  9.81  m/s2,  and  T  =  288.15°K. 

The  flat  plate  was  run  at  the  same  Reynolds’  numbers  as  the  airfoil.  Therefore, 
the  length  of  the  plate  will  be  set  to  0.202  meters  to  accommodate  a  similar  velocity, 
density,  and  viscosity  as  the  airfoil  experiment.  For  verification  that  the  viscosity  is  in 
the  correct  range,  the  values  for  air  at  STP  are  referenced  in  Table  1. 
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Table  1.  Standard  Temperature  and  Pressure  (STP)  for  Air 


Temperature: 

288.15  °K 

Pressure: 

101325  Pa 

Speed  of  Sound: 

340.2  m/s 

ft: 

1.7894x10  5N  •  s/m2 

p: 

In  White  [11],  it  is  found  that  these  flows  will  be  laminar,  not  turbulent,  in  nature 
when  comparing  the  Reynolds  numbers  to  the  statement  “the  boundary- layer  flow  is 
likely  to  be  laminar  in  the  range  1000  <  Re  <  106”  [1 1 :218]. 

For  the  two  flow  speeds,  the  maximum  and  minimum  boundary  layer  thickness 
was  used  to  facilitate  creation  of  a  flat  plate  grid.  The  value  of  x  was  set  to  75%  of  the 
length  of  the  plate,  x=0. 1515m,  to  ensure  that  the  boundary  layer  is  fully  developed  and 
to  simulate  where  the  DBD  will  be  placed  along  the  airfoil  as  seen  in  Figure  16.  For 
further  validation  and  calculations,  a  simulation  of  the  boundary  layer  velocity  profile 
was  performed  for  a  value  of  x  set  to  near  100%  of  the  length  of  the  plate,  x  =  0.2015  m. 
A  simulation  at  this  point  allowed  for  easy  verification  calculations  with  the  Fluent® 
simulations  and  for  determining  the  maximum  height  needed  for  the  grid  spacing. 

To  examine  how  the  boundary  layer  will  behave  along  the  flat  plate  at  different 
velocities,  the  following  equation  sets  in  Table  2  as  well  as  Figure  31  and  Figure  32  were 
created.  The  graphs  and  equation  sets  depict  a  boundary  layer  velocity  profile  at  a 
position  75%  and  near  100%  down  the  length  of  the  plate  as  already  described.  The 
analytic  equations  used  are  the  Blasius  boundary  layer  solutions  to  the  differential  flow 
equations. 
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Table  2.  Blasius  Boundary  Layer  Velocity  Profile  Calculations 

v  =  [/  =  2  m/s 


p  =  1 .225  kg/m3  (i  =1.7910“5N ■  s/m2  u  =  ^  =  1 .46 x  10“5 m2/s 

P 


=  (75%)  (0.202)  =0.15 15  meters 


Re  =  x  =  -x  =  0.021  xlO6 

\x  x> 


x  =  0.2015  meters 
Re  =0.027  xlO6 


'99 


5.5* 


=  5.79x10  meters  =  5.79  mm  8 


99 


5.5x 


=  6.67x10  ’meters  =  6.67  mm 


2  \ 

2 y _ 

^  Sv  2 

y  °99  °99  / 


=  691y-59700y': 


u~Un 


2  'l 

2y _ y_ 

^  Sv  2 

y°99  °99  J 


=  599y-  44900y2 


Figure  3 1 

2  m/s  Freestream  Velocity  Profile  at  x=0.1515  meters 
Analytic  Blasius  Velocity  Profile  VS  Fluent  Data 
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— * — Blasius  Profile  x=0.2015m  — ■ — Fluent  x=0.2015m 
Figure  32 

2  m/s  Free  stream  Velocity  Profile  at  x=0.2015  meters 
Analytic  Blasius  Velocity  Profile  VS  Fluent  Data 

Horizontal  Grid  Spacing 

The  flat  plate  grid  is  the  main  test  platform  and  was  modeled  prior  to  the  airfoil 
grid  creation.  The  flat  plate  used  a  total  of 250  cells  from  0  to  0.202  meters  with  a 
tanh(x)  spacing  and  initial  spacing  of  1  x  10  4  at  the  leading  edge.  A  total  of  96  points 
from  0  <  x  <  0.55c  (0  <  x  <  0.1 1 1)  meters  was  set;  where  the  cord  length  is  defined  as 

c=0.202m  from  the  table  in  Figure  16.  The  DPIV  data  depicted  in  Figure  3  was  used  as 
a  guide  for  increased  cell  density  as  it  shows  significant  wall  jet  induced  velocities  from 
(-0.020  <d  <  +0.060)  meters,  where  d  =  0.75c  ( d  =0.1515  m)  from  Figure  16. 
Finally,  small  cell  sizes  yielding  more  fidelity  but  more  calculation  time  and  possible 
instabilities  were  used  in  the  DBD  profile  region 
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The  largest  cell  in  the  region  of  the  DBD  depends  upon  the  force  density  scheme 
applied.  Thus,  a  grid  for  the  Boeuf  and  Pitchford  and  a  separate  grid  for  the  Roy  and 
Gaitonde  profiles  were  created.  The  Boeuf  and  Pitchford  profile  is  so  incredibly  small 
compared  to  the  entirety  of  the  grid  that  101  points  were  compressed  into  1mm  between 

(0.151m  <  x  <  0.152m) .  The  even  spacing  of  these  points  provides  a  distance  of  10  (im 
per  side  of  the  cell.  The  Roy  and  Gaitonde  profile  is  much  larger,  and  therefore  has  a 
span  of  39mm  between  (0.141m  <  x  <  0.180m) .  The  even  spacing  of  157  points  yields  a 
distance  of  250  (im  per  side  of  the  cell. 
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Vertical  Grid  Spacing 

The  maximum  height  of  the  flat  plate  grid  needs  to  be  calculated  using  the 
thickest  boundary  layer  expected,  which  is  obtained  when  the  free  stream  velocity  is 

slowest.  The  maximum  height  is  set  greater  than  or  equal  to  4  ■  (S,,^  ^  ) ,  where  899s| 

was  originally  set  for  a  15.2  m/s  scenario;  this  yielded  4-  (2.62  xlO~3)  meters,  or 

10.2mm.  The  height  of  the  smallest  cell  needs  to  be  calculated  in  just  the  opposite 
manner,  when  the  thinnest  boundary  layer  is  expected,  which  is  when  the  free  stream 

899 

velocity  is  fastest.  The  minimum  height  is  set  to  less  than  or  equal  to  — in  order  to 
have  about  10  cells  in  the  y-direction  to  capture  the  boundary  layer,  where  S99p  bl  was 


originally  set  for  a  30.4  m/s  scenario;  this  yielded 


1.60  xlO~3 
10 


meters,  or  0.160mm. 


Again,  force  density  profiles  play  a  role  in  the  size  of  the  vertical  grid  spacing. 
Boeuf  and  Pitchford  data  require  a  very  limited  area  for  their  force  density  profile,  and  as 
such  the  grid  to  be  used  for  simulating  their  profile  is  built  so  that  5 1  points  are  contained 
in  the  first  0.1mm,  a  spacing  of  2  |tm.  The  Roy  and  Gaitonde  profile  is  again,  much 
larger,  and  extends  to  1cm  off  of  the  surface  of  the  flat  plate.  A  total  of  101  points  are 
contained  in  this  region  for  a  spacing  of  100  (im  The  cell  spacing  for  either  of  these 
profiles  is  less  than  the  smallest  cell  size  as  determined  by  the  fastest  free  stream  velocity. 

Each  cell  height  above  the  profile  areas  follows  a  tanh(x)  spacing  approach  in 
Gridgen®.  The  smallest  cells  are  at  the  bottom  of  the  grid  where  the  DBD  is  simulated 
and  the  largest  cells  are  at  the  top  of  the  grid  where  the  free  stream  velocity  is. 
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0.102m  _ _ oin:'n 


The  flat  plate  grids  were  constructed  with  the  specifications  for  vertical  and 


horizontal  grid  point  spacing  as  seen  in  Figure  33  and  Figure  34. 


o  ;o:m 
Figure  33 

Boeuf  and  Pitchford  Flat  Plate  Grid  Geometry 


0.202m 

Figure  34 

Roy  and  Gaitonde  flat  Plate  Grid  Geometry 
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Flat  Plate  Validation 


The  Flat  Plate  grids  for  the  Boeuf  and  Pitchford  profile  and  the  Roy  and  Gaitonde 
profile  were  validated  by  examining  the  boundary  layer  profile  against  a  Blasius  profile. 
To  obtain  the  boundary  layer  velocity  profile  for  each  of  the  grids,  an  unsteady  viscous 
simulation  was  run  using  0.001  second  time  steps  for  1000  steps  (1  second  total  time)  at 
2  m/s  ,  15.2  m/s ,  and  30.4  m/s .  For  each  time  step,  a  maximum  of  20  sub -iterations 
could  be  performed  before  moving  to  the  next  time  step.  The  residual  was  also 
monitored  for  convergence  to  10‘6  for  each  time  step,  a  condition  that  would  cause  a 
move  to  the  next  time  step  prior  to  20  sub -iterations  being  completed.  Further  details 
concerning  the  simulation  setup  can  be  found  in  Appendix  C. 
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NACA  0009  Airfoil  Grid  Geometry 


The  creation  of  the  NACA  0009  Airfoil  grid  involved  three  data  sets,  each 
increasing  in  fidelity.  The  grids  were  constructed  in  Gridgen®.  Airfoil  Data  Set  #1  and 
#2  found  in  Table  11  and  Table  12,  respectively,  in  Appendix  A,  were  first  used  in  an 
attempt  at  validation  However,  neither  was  found  to  have  the  fidelity  needed  for  the 
NACA0009  airfoil  shape  at  the  Leading  Edge  (LE).  The  next  method  explored  is  a 
mathematical  formula  [13]  which  yields  much  higher  fidelity  in  this  region 


Given  a  4-digit  symmetric  airfoil  such  as  the  NACA  0009,  we  can  decipher  its 


naming  convention  to  give: 


NACA  0009 
NACA  f  xf  t  t 


f  maximum  camber 
Xf  position  of  maximum  camber 
x  position  along  x-axis 

t  choriiength  >  the  digits  represent  a  %,  therefore  09=9%=.09 


To  locate  a  2-D  coordinate,  the  following  equations  are  employed. 


^  x  +  yt  •  sinO  ^ 

yc  ±  yt ■ cos0 


where 


(30) 


tanO  = 


dyc 

dx 


(31) 


5T  =  f _ L_ 

c  c(l-^)2 


(l-  2x:)  +  2x, : — 
c 


Vc  J 


(32) 


—  =  5t  [0.29690x° 5  -  0. 12600x  -  0.35 160x2  +  0.28430x3  -  0. 10150x4]  (33) 
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With  no  camber,  f=0  and  Xf=0,  the  above  equations  simplify  to: 


where 


f X ) 

(  x  ) 

+i 

(34) 


yt=c-5t[0.29690/'5- 0.1 2600x- 0.35 160.r  +  0.28430.V3 -0.10150*4]  (35) 

These  equations  will  be  used  in  the  UDF  to  find  points  along  the  airfoil,  as  well  as  give 
the  tangent  slope  at  a  given  x- point  by  deriving  the  following  equation  from  the  above. 

yt'  =  c  -  5r[0.5*0.29690*~°5  -0.12600-2*0.35160x+3*0.28430x2  -4*0.10150x3  J  (36) 

The  resulting  airfoil  shape  is  much  smoother  than  the  previous  data  sets  given. 

The  data  points  are  listed  in  Appendix  A.  At  the  leading  edge  (LE),  the  difference  is  the 
most  dramatic.  At  the  trailing  edge  (TE),  the  difference  is  negligible.  The  overall  shape 
of  the  airfoil  has  not  changed  by  employing  an  equation  for  a  model.  However,  more 
fidelity  was  gained  in  specific  regions,  such  as  the  leading  edge.  A  visual  comparison 
between  the  data  sets  at  the  leading  and  trailing  edges  was  performed  in  Figure  35  and 
Figure  36.  The  increase  in  data  points  has  given  rise  to  a  much  smoother  curve  in  both  of 
these  areas. 


in-ii 


meters 


— ♦ —  Equation  M  Airfoil  Data  Set  #1  ^  Airfoil  Data  Set  #2 

Figure  35 

Leading  Edge  of  Airfoil  Comparison  between  Equation  and  Data  Sets 


meters 


— ♦ —  Equation  U  Airfoil  Data  Set  #1  —A—  Airfoil  Data  Set  #2 

Figure  36 

Trailing  Edge  of  Airfoil  Comparison  between  Equation  and  Data  Sets 


The  airfoil  grid  was  constructed  using  the  parameters  of  grid  spacing  for  the  Roy 
and  Gaitonde  flat  plate  profile.  The  top  section  was  first  constructed  and  then  mirrored  to 
create  a  symmetric  bottom  section.  The  top  section  was  then  increased  in  its  number  of 
points  and  its  horizontal  grid  spacing  was  altered  to  fit  the  horizontal  grid  spacing  of  the 
flat  plate.  The  bottom  section  of  the  airfoil  remained  sparser  in  points  as  there  was  no 
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need  for  a  greater  fidelity  along  the  cord  of  the  airfoil.  Points  were  clustered  via  a  “tanh” 
function  along  the  airfoil  curve  towards  the  leading  edge  for  both  sections  of  the  airfoil, 
where  points  steadily  increase  in  their  separation  distance  for  smooth  cell  size  changes. 
Finally,  an  elliptical  solver  was  run  in  Gridgen®  to  make  the  cells  closest  to  the 
boundaries  orthogonal.  The  “tanh”  spacing  and  the  elliptical  solver  have  the  effect  of 
aiding  in  solution  stability. 

The  final  grid  geometry  is  displayed  in  Figure  37,  showing  the  increased  density 
of  points  near  the  boundary  layer  and  on  the  top  of  the  airfoil  at  75%  cord  length  where 
the  DBD  is  to  be  simulated. 


Figure  37 

Final  NACA  0009  Airfoil  Grid  Geometry 
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NACA  0009  Airfoil  Validation  Confirmation 


Experimental  data  was  obtained  fromSelig  [15]  and  is  compared  to  simulated 
data  obtained  for  validation  of  the  grid  mesh  Data  labeled  SATurb  or  Laminar  preceding 
the  Reynold’s  number  is  Fluent®  simulated  data.  A  Spalart-Allmaras  Turbulence  model 
was  employed  to  obtain  the  data  shown  in  Figure  38  thru  Figure  41  and  in  additional  data 
presented  in  Appendix  A  in  Figure  69  thru  Figure  72.  As  seen  in  Figure  38  thru  Figure 
41,  the  Laminar  model  data  curves  did  not  match  as  well  as  the  turbulence  model  data 
when  compared  to  experiments. 
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AoA  vs  CL 


-b— Re=1 50600  -e— Re=1 99900  Laminar  Re=1 80000  SATurb  Re=1 80000 

Figure  38 
Rec=l  80,000 

NACA  0009  Angle  of  Attack  VS  Lift  Coefficient 
Comparison  of  Laminar  Model  and  Spalart-Allmaras  Turbulence  Model 
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— • — Laminar  Re=1 80000  — * — SATurb  Re=1 80000 


Figure  39 
Rec=l  80,000 

NACA  0009  Angle  of  Attack  VS  Lift  Coefficient  Difference 
Comparison  of  Laminar  Model  and  Spalart-Allmaras  Turbulence  Model 
as  Compared  to  Experimental  Data  by  Selig  [15] 
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AoA  vs  CL 


— b—  Re=302900  —•—Laminar  Re=360000  — *—  SATurb  Re=360000 

Figure  40 
Rec=360,000 

NACA  0009  Angle  of  Attack  VS  Lift  Coefficient 
Comparison  of  Laminar  Model  and  Spalart-Allmaras  Turbulence  Model 


AoA  vs  CL 
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— • — Laminar  Re=360000  — •—  SATurb  Re=360000 


Figure  41 
Rec=360,000 

NACA  0009  Angle  of  Attack  VS  Lift  Coefficient  Difference 
Comparison  of  Laminar  Model  and  Spalart-Allmaras  Turbulence  Model 
as  Compared  to  Experimental  Data  by  Selig  [15] 
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The  physical  data  that  is  displayed  in  the  previous  figures  (Figure  38  thru  Figure 
41)  was  not  tripped  to  turbulent  over  the  airfoil  for  the  purpose  of  the  measurements  [15]. 
This  was  confirmed  with  Selig  [15]  via  e-mail  correspondence.  As  a  result,  it  is  expected 
that  a  Laminar  model  should  be  applied.  Further,  turbulence  is  usually  calculated  to 

occur  above  a  Reynolds’  number  of  0.5  xlO6 .  However,  there  is  indication  that 
turbulence  and  separation  are  occurring  at  the  half-cord  on  this  symmetric  airfoil.  The 
Reynolds’  numbers  may  be  approaching  a  lower  limit  for  consideration  of  a  turbulence 
model  such  as  Spalart-Allmaras.  Figure  39  and  Figure  41  show  that  the  Spalart-Allmaras 
data  is  much  more  consistent  with  Selig’ s  experimental  data  as  compared  with  the 
laminar  data. 

The  turbulent  simulated  data  obtained  more  closely  mirrors  the  physical  data 
given  by  independent  sources  [15].  This  suggests  that  the  model  is  best  represented  using 
the  Spalart-Allmaras  Turbulence  Model  in  Fluent®.  However,  because  turbulence  was 
not  tripped  during  the  actual  physical  measurements,  this  validation  can  not  be 
substantiated  and  a  laminar  solver  is  expected  to  be  run.  Nonetheless,  the  laminar  solver 
will  not  be  run  due  to  the  poor  behavior  of  the  model  as  compared  to  measured 
experimental  data.  Instead,  for  purposes  of  this  simulation,  it  is  assumed  that  the  Spalart- 
Allmaras  turbulence  model  will  apply  so  as  to  give  data  that  more  closely  mirrors  reality. 
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IV.  Results 


It  was  necessary  to  examine  different  pieces  of  the  DBD  simulation’s  behavior  in 
order  to  arrive  at  a  parameterized  view  of  the  DBD’s  operation.  Several  data  sets  were 
simulated  and  results  returned.  Thermal  energy  addition  versus  momentum  addition 
effects  were  simulated  as  well  as  power  per  unit  length  effects.  This  section  will  cover 
these  simulation  sets  and  their  results. 

Processing 

The  airfoil  data  sets  simulated  required  a  great  deal  of  processing  power.  Each  set 
of  data  for  the  airfoil  was  simulated  at  angles  of  attack  from  0°  to  +16°  in  2°  increments 
for  each  of  the  2  Reynolds’  numbers  used  in  validation;  a  total  of  18  simulations  for  each 
set.  The  simulations  were  processed  on  the  AFIT  computer  cluster  machine  Tahoe. 

Tahoe  is  a  64  node  computer  cluster  with  each  node  having  2  AMD  Opteron 
2.2GHz  processors.  The  node  set  of  processors  reside  on  the  same  control  board  and 
have  4GB  of  RAM  shared  in  a  NUMA,  Non-Uniform  Memory  Architecture,  where  each 
processor  on  that  node  is  guaranteed  its  part  of  the  RAM,  not  shared.  Connections 
between  the  nodes  are  handled  using  lGbps  Ethernet. 

AFIT  has  ownership  of  30  Main/Startup  Fluent®  licenses  with  88  Multi-processor 
Fluent®  licenses.  The  main  Fluent®  licenses  along  with  available  nodes  drove  the 
overall  rate  at  which  each  set  was  completed.  Each  simulation  was  performed  on  1  node, 
using  both  processors  in  an  attempt  to  reduce  computational  time.  The  average 
completion  time  for  a  simulation  was  12  hours,  with  each  set  taking  24  hours  to 
complete.  The  validation  cases  and  another  set  were  run  using  only  single  processors, 
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which  were  found  to  be  1.5  times  faster.  This  was  discovered  after  most  of  the 
simulations  had  been  run  and  is  not  a  surprising  result  for  such  a  small  2D  simulation. 
The  decreased  performance  with  increased  processor  numbers  is  due  to  the  network 
transfer  speed  being  a  significant  factor  in  the  processing  of  each  iteration;  which  means 
the  processors  are  idle  while  they  transfer  information  between  each  other  and  are 
therefore  less  efficient.  Future  work  sho  uld  use  only  a  single  processor  for  each 
simulation  as  a  result. 

Appendix  C  details  the  PBS  (Portable  Batch  System)  script  and  journal  files  that 
allowed  automation  of  each  simulation. 

Matching  Force  Density  Profiles 

Several  trials  were  conducted  in  order  to  fit  the  force  density  profiles  similar  to 
those  reported  by  either  Boeuf  and  Pitchford  in  Figure  18  and  Figure  19  or  Roy  and 
Gaitonde  in  Figure  23  and  Figure  24.  This  was  due  to  the  boundary  layer  velocities 
affecting  the  profiles  after  steady  state  had  been  reached.  As  was  discussed  in  the 
Simulation  Setup  section,  an  increase  in  the  force  density  increases  the  velocity  of  the 
flow,  but  the  velocity  increase  is  indirectly  proportional  to  the  flow  velocity  where  the 
force  density  is  acting.  The  boundary  layer  has  an  increasing  velocity  as  the  distance 
from  the  wall  is  increased;  corresponding  to  a  force  density  that  drops  as  the  distance 
from  the  wall  increases. 

Velocity  profile  data  reported  from  these  simulations  gives  its  location  relative  to 
the  overlap  between  the  electrode  transition,  also  reported  as  the  center  of  the  DBD. 
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Boeuf  and  Pitchford  Force  Density  Profile 

The  Boeuf  and  Pitchford  profile  has  no  appreciable  effects  on  the  flow  at  low 
velocities  of  2  m/s,  and  as  a  result,  less  of  an  effect  at  higher  velocity  flows.  Figure  43  to 
Figure  46  show  that  within  the  boundary  layer  there  is  minimal  deviation  from  the 
baseline  boundary  layer  velocity  profile.  A  wall  jet  resulting  in  deviation  from  the 
baseline  case  is  expected  as  far  away  as  75mm  from  the  DBD  source,  as  is  seen  in  Figure 
3.  The  electrode  length  in  Figure  3  is  2  orders  of  magnitude  greater  than  that  of  the 
Boeuf  and  Pitchford  simulation.  The  simulated  DBD  force  density,  shown  in  Figure  42, 
must  be  greater  in  this  spatially  confined  profile  in  order  to  be  effective  in  its  confined 
volume. 


Table  3.  Boeuf  and  Pitchford  Matching  Data  Set 


Power  (W) 

1  meter 

Thermal 
Energy  % 

Thermal 
Power  3L 

Momentum 

Power 

X- Momentum 

Y-M 

omentum 

% 

Power  ^ 

% 

Power 

0.0003 

10% 

0.00003 

0.00027 

100% 

0.00027 

0% 

0 

uam ■?.  SOO  1000  1500  3000  3S00  3000  3500  4000  <500  5000  5500  6000  6500  7000  7500  9000  BBOO  9000  9500  10000  10500 


Figure  42 

X- Momentum  Force  Density  for  the  Case  in  Table  3 
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Figure  43 

Boeuf  and  Pitchford  Weighting  Profile  2  m/s  Simulation  Result  at  x=0mm 

Compared  to  Baseline 
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Figure  44 

Boeuf  and  Pitchford  Weighting  Profile  2  m/s  Simulation  Result  at  x=0mm 
Compared  to  Baseline  (Close-up) 
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Figure  45 

Boeuf  and  Pitchford  Weighting  Profile  2  m/s  Simulation  Result  at  x=5mm 

Compared  to  Baseline 
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Figure  46 

Boeuf  and  Pitchford  Weighting  Profile  2  m/s  Simulation  Result  at  x=5mm 
Compared  to  Baseline  (Close-up) 
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Several  cases  were  examined  using  the  force  density  profile  of  an  extremely 
confined  volume  by  which  to  represent  the  DBD.  Figure  48  displays  these  simulation 
results  and  compares  against  the  baseline  case  at  x=5mm  from  the  DBD  location.  The 
simulations’  specific  settings  are  listed  in  Table  4.  These  settings  resulted  in  the  force 
density  shown  in  Figure  47  and  represent  a  2  order  of  magnitude  increase  from  the  values 
suggested  by  Boeuf  and  Pitchford. 


Table  4.  Boeuf  and  Pitchford  Data  Set  for  Increasing  Thermal  Addition  with 
Fixed  Momentum  Percent  Addition 


Power  (W) 

1  meter 

Thermal 
Energy  % 

Thermal 
Power  ^ 

Momentum 
Power  ^ 

X-Momentum 

Y-M 

omentum 

% 

Power 

% 

Power 

5.0 

40% 

2.0 

3.0 

70% 

2.10 

30% 

0.90 

50% 

2.5 

2.5 

1.75 

0.75 

60% 

3.0 

2.0 

1.40 

0.60 

70% 

3.5 

1.5 

1.05 

0.45 

Figure  47 

X-Momentum  Force  Density  for  40%  Thermal  Energy  %  in  Table  4 
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Figure  48 

Boeuf  and  Pitchford  Force  Density  Simulation  Result  at  x=5mm 
Compared  to  Baseline  for  Test  Cases  in  Table  4 
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As  the  amount  of  thermal  energy  was  increased  and  therefore  the  momentum 
energy  decreased,  the  boundary  layer  profiles  in  Figure  48  decrease  towards  the  baseline 
case  as  expected.  The  maximum  velocity  magnitude  is  found  at  the  DBD  source  for 
these  simulations  reaching  close  to  6  m/s,  but  rapidly  falls  back  to  the  baseline  velocity 
profile.  What  was  expected  from  Boeuf  and  Pitchford’s  original  force  density  profile  was 
a  flow  with  a  wall  jet  that  has  a  maximum  velocity  of  approximately  8  m/s  to  10  m/s  as 
was  shown  in  Figure  20  or  a  minimum  of  2.3  m/s  as  shown  in  Figure  28.  Boeuf  and 
Pitchford  [16]  state  that  their  time  integrated  force  density  profile  will  “...approximately 
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give  the  contours  of  constant  velocity  increment  in  units  of  m/s”.  There  is  a  discrepancy 
between  the  ~8  m/s  wall  jet  predicted  by  Boeuf  and  Pitchford,  the  1.9  m/s  predicted  with 
stationary  flow  energy  addition  analysis  in  the  Simulation  Setup,  and  the  0.15  m/s  that 
occurred  in  the  Fluent®  simulation  The  discrepancy  cannot  be  explained  at  this  time. 

In  a  final  attempt  to  explain  this  velocity  discrepancy,  a  final  set  of  simulations 
were  performed  to  examine  the  added  effect  of  increasing  the  extent  of  the  force  density 
profile  in  the  horizontal  direction.  The  four  cases  examined  are  listed  in  Table  5.  It  was 
hypothesized  by  Boeuf  and  Pitchford  [16]  that  their  simulation  space  was  not  large 
enough  and  that  by  increasing  its  horizontal  extent,  the  force  density  would  continue. 

This  is  a  conclusion  that  can  be  drawn  from  examining  both  impulse  density  results  for 
their  400pm  and  800pm  wide  DBD  simulations;  they  behave  similarly  but  with  a  larger 
horizontal  extent  for  the  force  density  when  examining  the  800pm  case.  Therefore,  if  the 
horizontal  extent  of  the  DBD’s  force  density  was  allowed  to  increase,  the  flow  would 
undergo  more  accumulated  impulse  and  have  a  larger  velocity  increase.  Initially,  in 
Figure  50,  the  velocities  begin  with  near  the  same  wall  jet  velocity  increase  at  x=0mm. 

In  the  following  figures.  Figure  51,  Figure  52,  and  Figure  53,  velocity  profiles  are  taken 
an  additional  5mm  from  x=0mm  for  each  subsequent  figure.  It  is  apparent  that  there  is  a 
small  wall  jet  velocity  increase  due  to  the  lengthened  horizontal  spatial  extent.  However, 
as  soon  as  the  velocity  profile  is  no  longer  being  measured  within  a  force  density  profile, 
the  wall  jet  abruptly  returns  to  baseline  velocity  profile  behavior.  This  is  still  not  the 
suggested  behavior.  It  may  be  possible  that  the  vertical  extent  of  the  force  density  will 
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also  need  to  be  increased,  such  that  a  tendency  is  more  towards  a  force  density  profile  of 


Roy  and  Gaitonde’s  reports. 


Table  5.  Boeuf  and  Pitchford  Data  Set  for  Increasing  Simulation  Extent  while 
maintaining  Force  Density  Profiles  Suggested  by  Boeuf  and  Pitchford 


Power  (W) 

1  meter 

Simulation 
Extent  (mm) 

0.0003 

0.8 

0.0005 

1 

0.0025 

5 

0.0050 

10 

-&■ 

■& 

-<s- 

■v- 


Baselme 

0.0003W 

0.8mm 

0.0005W 

limn 

0.0025W 

5imn 

0.0050W 

lOirun 

0  4  o- Thermal 
0°  o- Thermal 
0®/o- Thermal 
0°  o-Thennal 


100*  o-Momentum 
100®  o-Momentuin 
100®  o-Momentiun 
100°  o-Momentiun 


100®  «-X  Momentiun 
100®  «-X  Momentiun 
100®  o-X  Momentiun 
100®  «-X  Momentiun 


0®  b-Y  Momentiun 
0®  o-Y  Momentiun 
0%-Y  Momentiun 
0®  o-Y  Momentiun 


Figure  49 

Fegend  for  the  Velocity  Profiles  Fisted  in  Table  5  and 
Displayed  in  Figure  50  to  Figure  53  (Close-up) 


Boeuf  and  Pitchford  2m/s  Increasing  Simulation  Extent  Velocity  Profiles  at  x=0mm 
Compared  to  Baseline  for  Test  Cases  in  Table  5  (Close-up) 
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Boeuf  and  Pitchford  2m/s  Increasing  Simulation  Extent  Velocity  Profiles  at  x=5mm 
Compared  to  Baseline  for  Test  Cases  in  Table  5  (Close-up) 


Boeuf  and  Pitchford  2m/s  Increasing  Simulation  Extent  Velocity  Profiles  at  x=10mm 
Compared  to  Baseline  for  Test  Cases  in  Table  5  (Close-up) 
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Boeuf  and  Pitchford  2m/s  Increasing  Simulation  Extent  Velocity  Profiles  at  x=15mm 
Compared  to  Baseline  for  Test  Cases  in  Table  5  (Close-up) 


Roy  and  Gaitonde  Force  Density  Profile 

The  Roy  and  Gaitonde  profile  has  effects  at  a  free  stream  velocity  of  2  m/s.  This 
section  details  the  required  settings  and  results  of  the  Roy  and  Gaitonde  force  density 
profile  on  a  flat  plate  with  a  free  stream  velocity  of  2  m/s. 

The  setup  that  allows  for  the  given  profile  is  listed  in  Table  6.  The  two  figures 
that  follow  the  table,  Figure  54  and  Figure  55,  show  that  the  force  density  profile  was 
simulated  correctly  with  the  settings  in  Table  6.  This  is  not  entirely  evident  in  Figure  54 
due  to  the  negative  region  at  0.16m  having  less  of  a  peak  absolute  magnitude,  and 
therefore  not  showing  in  the  same  fidelity  as  the  positive  regions. 

The  wall  jet  is  evident  at  x=5mm  as  seen  in  Figure  56  and  continues  downstream 
in  an  accelerating  manner  to  a  peak  of  4.2  m/s  at  10mm  from  the  DBD  source.  The  wall 
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jet’s  velocity  magnitude  then  starts  to  fall  off  due  to  fluid  sheer,  but  remains  near  4  m/s  at 
50mm  from  the  DBD  source.  This  behavior  is  nearly  twice  the  expected  velocity 
estimated  in  the  Simulation  Setup  section  as  well  as  depicted  in  Figure  25  with 
experimental  DBD  setup  and  measurements,  but  fulfills  the  original  goal  of  producing  a 
near  5  m/s  wall  jet. 


Table  6.  Settings  to  Achieve  Roy  and  Gaitonde  Weighting  Profile 


Power  (W) 

1  meter 

Thermal 
Energy  % 

Thermal 
Power  2L 

Momentum 
Power  ^ 

X- Momentum 

Y-M 

omentum 

% 

Power  ^ 

% 

Power 

5.0 

97.8% 

4.89 

0.11 

90% 

0.099 

10% 

0.011 

N_ 

m3 


udm  ?.  -1000  -500 


500  IOCO  1500  7000  7500  3000  3500  4000 


Figure  54 

Roy  and  Gaitonde  X- Momentum  Force  Density  Profile  from  Table  6 
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Figure  55 

Roy  and  Gaitonde  Y-Momentum  Force  Density  Profile  from  Table  6 


■B -  5W  97.8%- Thermal  2.2%-Momentum:  90%-X Momentum  10 %-Y Momentum 

Figure  56 

Boundary  Layer  Velocity  Magnitude  Profile  for 
Roy  and  Gaitonde  Force  Density  Simulation  Result  at  x=5mm 
Compared  to  Baseline  for  Test  Case  in  Table  6 
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Thermal  Energy  Dependence 

It  has  been  noted  by  Newcamp  [18],  that  there  is  a  limit  to  effectiveness  of  the 
DBD  input  power.  Newcamp  found  that  as  the  power  was  increased,  and  the  flow 
velocity  remained  unchanged,  the  velocity  of  the  wall  jet  decreased.  The  suspect 
mechanism  is  that  the  amount  of  momentum  transfer  to  the  flow  remains  unchanged 
while  the  thermal  energy  addition  increases.  For  a  gaseous  fluid,  increasing  its 
temperature  will  result  in  a  fluid  that  is  more  viscous,  or  more  difficult  to  move.  This  is 
the  opposite  effect  that  temperature  increase  will  have  on  liquid  fluids. 

To  confirm  or  deny  the  effects  of  thermal  energy  being  the  mechanism  for 
slowing  down  the  jet,  several  simulations  were  run  that  fixed  the  momentum  sources  in 
both  the  x  and  y  directions,  yet  allowed  for  the  thermal  energy  addition  to  be  varied  in  a 
2  m/s  freestream  flow  over  a  flat  plate.  The  specific  settings  are  displayed  in  Table  7. 
The  Roy  and  Gaitonde  force  density  profile  was  used,  as  it  has  appreciable  effects  on  the 
boundary  layer  flow  as  compared  to  the  Boeuf  and  Pitchford  force  density  profile. 

The  high  and  the  low  end  of  input  power  for  Table  7’s  range  was  derived  using 
the  near  minimum  amount  of  power  it  would  take  to  create  an  approximate  5  m/s  flow 
and  the  maximum  that  was  put  into  a  DBD  before  it  physically  failed  and  could  never  be 
used  again  from  Newcamp’s  experiments  [18].  Newcamp  reports  that  failure  occurred  at 
25W/5"  =  200  W/m . 

The  boundary  layer  ends  6mm  above  the  surface  of  the  flat  plate  for  a  2  m/s  flow. 
The  boundary  layer  profile  and  the  profiles  of  the  simulations  run,  listed  in  Table  7,  are 


IV- 14 


displayed  in  Figure  57.  As  expected,  Figure  57  shows  that  as  the  thermal  energy  is 
increased,  the  wall  jet  velocity  will  decrease.  However,  it  is  also  seen  that  this  decrease 
is  negligible,  varying  velocity  from  4.6  m/s  to  4.4  m/s,  a  0.2  m/s  difference  that  decreases 
as  the  thermal  power  addition  is  increased  over  4  orders  of  magnitude.  Further,  Figure  59 
shows  that  the  linear  temperature  increase  seen  is  consistent  with  estimated  calculations 

Ae  _2, 

using  the  equation  AT  = - ;  where  kb  - 1 .380658  10  is  the  Boltzmann  Constant, 

jkbN 

Ae  =  Power  -a t  where  a t  is  the  user  defined  time  step  taken  by  Fluent®  to  resolve  a  time 

,  .  ,  A7  1  mole  6.0221367  TO23#  , 

accurate  solution,  and  N  =  p - vol . 

0.029kg  1  mole 

This  simulation  set  from  Table  7  does  account  for  a  decrease  in  velocity 
magnitude  as  the  input  power  is  increased;  as  described  by  Newcamp  [18].  However,  to 
be  consistent  with  Newcamp’ s  reporting  in  Figure  58,  a  more  accurate  scale  to  compare 
on  would  have  Figure  57’s  data  divided  by  U„,  the  fluid’s  free  stream  velocity.  As 
=  2  m/s  ,  the  peak  magnitudes  would  then  read  2.3  and  2.2  for  the  lower  and  higher 
thermal  energy  cases,  respectively.  This  is  a  0.1  difference  that  is  similar  to  the  data 
presented  in  Figure  58. 
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Table  7.  Thermal  Energy  Dependence  Test  Case  Sets  for  Roy  and  Gaitonde  Profile 


Power  (W) 

Thermal 

Thermal 

Momentum 

X- Momentum 

Y- Momentum 

Power 

w 

m 

1  meter 

Energy  % 

Power  ^ 

Power  ^ 

% 

Power  ^ 

% 

0.2 

10% 

0.02 

1 

82% 

0.82 

2 

91% 

1.82 

3 

94% 

2.82 

4 

95.5% 

3.82 

5 

96.4% 

4.82 

10 

98.2% 

9.82 

0.18 

55% 

0.099 

45% 

0.011 

15 

98.8% 

14.82 

20 

99.1% 

19.82 

25 

99.28% 

24.82 

80 

99.775% 

79.82 

120 

99.85% 

119.82 

160 

99.8875% 

159.82 

200 

99.91% 

199.82 

0.006 


0.004 


>- 

0.002 


velocity-magnitude  (m/s) 

-0 -  Baseline 

-A -  0.2W  10°  o- Thermal  90%-Moinentiun:  55 %-X  Momentum  45%-Y  Momentum 

-V -  200W  99.91®  o- Thermal  0.09%-Momentum:  55  %-X  Momentum  45%-Y  Momentmn 

Figure  57 

Velocity  Profile  on  Flate  Plate  5mm  Downstream  of  DBD  Upper  Electrode; 

Fixed  Momentum,  Varying  Thermal,  Rex=20.8k 


IV- 16 


Velocity  Profile  for  Re  =  10,000  at  7.1  mm  Aft  ofTop  Electrode 


DPIV  Velocity  Profile  7.1mm  Downstream  of  DBD  Upper  Electrode  [18]; 
Varying  Total  Power,  Rex=10k 


— ♦ —  Fluent  Simulated  Increase  — ■ — Estimated  Increase 
Figure  59 

Temperature  Increase  vs  Thermal  Power  Addition  for  Table  7  Simulation  Set 
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To  be  complete,  a  simulation  set  with  a  fixed  percentage  of  power  transferred  into 
thermal  addition  as  well  as  momentum  addition  with  varying  power  levels  was  run.  The 
simulation  set  from  Table  8  was  performed  on  the  basis  that  the  power  level  is  increased 
by  raising  the  voltage  of  the  system.  The  increased  voltage  does  cause  more  ionization  of 
the  fluid,  air,  and  therefore  a  higher  current  level  of  the  system.  Since  P  =  I  V  ,  this 
leads  to  the  increase  in  power  seen  in  the  experiments.  Because  the  current  is  increasing, 
it  is  not  unreasonable  to  assume  that  there  will  be  more  ionization  resulting  in  a  larger 
wall  jet  velocity  magnitude.  This  theory  would  result  in  the  opposite  of  what  is  seen  in 
Figure  58.  The  hypothesis  is  correct  in  stating  that  as  power  is  increased,  the  wall  jet 
velocity  magnitude  will  increase  as  is  shown  in  Figure  60.  The  largest  wall  jet  velocity 
magnitude  coincides  with  the  most  power  input  to  the  system,  with  the  smallest  wall  jet 
velocity  magnitude  coinciding  with  the  least  amount  of  power  input  to  the  system, 
respectively.  However,  this  is  not  what  is  seen  in  Figure  58,  and  as  such  this  partition 
must  be  discarded.  As  was  shown  in  the  thermal  energy  dependence  simulations  just 
prior  to  this  subsection,  there  is  minimal  effect  to  the  wall  jet  velocity  by  increasing  the 
thermal  energy  source  input  to  the  system.  However,  a  small  variance  in  the  momentum 
source  will  lead  to  a  significant  change  in  the  wall  jet  velocity. 


Table  8.  Roy  and  Gaitonde  Weighting  for  Varying  Power  Levels  with 
Constant  Percent  Momentum  and  Thermal  Addition 


Power  (W) 

1  meter 

Thermal 

Energy 

% 

Thermal 

Power 

w. 

m 

Momentum 

% 

Momentum 

Power 

w. 

m 

X-Momentum 

Y -Momentum 

% 

Power 

w 

m 

% 

Power 

w 

m 

80 

97.8% 

78.24 

2.2% 

1.76 

90% 

1.584 

10% 

0.176 

120 

117.36 

2.64 

2.376 

0.264 

160 

156.48 

3.52 

3.168 

0.352 

200 

195.60 

4.40 

3.960 

0.440 
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80W  97.8°  o-Tlieimnl  2.2%-Momentum:  90°  o-X  Moinentiun  10°  o-Y  Momentum 
120W  9". 8°  o-Tlieiinal  2.2%-Momentum:  90°  o-X  Moinentiun  10°  o-Y Momentum 


160W  9 
200W  9 


'.8°  «-Tliennal  2.2%-Momentum:  90°  «-X Moinentiun  10° o-Y Momentum 


8%- Thermal  2.2%-Momentum:  90 %-X Moinentiun  10 %-YT  Moinentiun 

Figure  60 

Boundary  Layer  Velocity  Magnitude  Profile  for 
Roy  and  Gaitonde  Weighting  at  x=5mm  for  Table  8 


Airfoil  Results 

The  airfoil  test  is  the  final  set  of  cases.  The  flat  plate  scenarios  have  allowed  for 
the  area  under  examination  to  be  narrowed  to  a  specific  weighting  profile  that  best 
represents  the  physical  effect,  as  well  as  how  to  setup  that  profile.  This  set  of  cases  will 
allow  for  insight  into  the  effects  a  DBD  may  have  upon  a  NACA0009  airfoil.  The  goal  is 
to  model  the  data  obtained  by  Corke  [7]  in  Figure  12,  Figure  13,  Figure  14,  and  Figure 
15.  The  force  density  profile  chosen  was  the  Roy  and  Gaitonde  profile  and  its  setup 
parameters  are  given  in  Table  9. 
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"able  9.  Roy  and  Gaitonde  Weighting  Profile 

Parameters  for  Airfoil 

Power  (W) 

1  meter 

Thermal 

Energy 

% 

Thermal 

Power 

_w_ 

m 

Momentum 

% 

Momentum 

Power 

w 

m 

X-Momentum 

Y -Momentum 

% 

Power 

w 

m 

% 

Power 

w 

m 

5.0 

97.8% 

4.89 

2.2% 

0.11 

90% 

0.099 

10% 

0.011 

The  airfoil  was  then  simulated  over  1000  time  steps  of  0.001  seconds  for  a  total  of 
1  second  simulation  time.  Up  to  20  iterative  steps  were  allowed  for  each  time  step  for 
increased  accuracy,  as  was  done  with  the  flat  plate  simulations .  The  setup  of  the  flow 
was  such  that  the  free  stream  velocities  were  15.2  m/s  and  30.4  m/s  to  yield  Re=180k  and 
Re=360k  at  0.75c  respectively. 

As  was  deduced  from  the  Corke  paper  [7],  the  placement  of  the  DBD  at  0.75c  has 
not  had  any  appreciable  effect  on  the  stall  characteristic  of  the  airfoil.  Similar  to  Figure 
12  and  Figure  14,  the  two  figures  generated  by  the  simulation  sets,  Figure  61  and  Figure 
62,  show  a  similar  minimal  effect. 


Figure  61 

NACA0009  Airfoil  Ci  vs.  AoA  for  Re=180k 


IV-  20 


The  opposite  is  the  case  when  comparing  Corke’s  results  [7]  of  Figure  13  and 
Figure  15  with  Figure  63  and  Figure  64  respectively.  The  Q  vs.  Cd  curves  are  not  as 
separated  and  distinct.  A  closer  examination  of  the  data  was  performed  by  taking  the 
difference  between  the  simulated  data  with  a  DBD  in  operation  and  the  baseline  cases 
without  a  DBD  in  operation. 


-0.1  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0 

Coefficient  of  Lift 


I — * — Baseline  SATurb  Re=180k  —  DBD  SATurb  Re=180k 

Figure  63 

NACA0009  Airfoil  Q  vs.  Cd  for  Re=180k 
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The  trend  of  a  positive  shift  is  still  seen,  but  it  is  very  small  as  shown  in  Figure  65 
and  Figure  66.  The  erratic  behavior  of  the  data  difference  above  12°  AoA  shows  the 
simulation  of  the  DBD  is  most  likely  only  valid  until  this  point.  Separation  of  the  flow  is 
suspected  to  have  occurred  at  this  angle  of  attack.  Because  the  DBD  is  located  at  0.75c, 
the  contribution  to  the  airfoil’s  performance  is  expected  to  be  minimal.  This  is  due  to  a 
minimal  amount  of  separation  occurring  at  about  0.5c  at  an  angle  of  attack  of  0°. 
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I  Re=1 80k  Difference  Re=360k  Difference! 

Figure  65 

NACA0009  Airfoil  Ci  Difference  for  DBD  [On-Off]  vs.  AoA  for  Re=180k 


Figure  66 

NACA0009  Airfoil  Q  Difference  for  DBD  [On-Off]  vs.  AoA  for  Re=360k 
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V.  Conclusions 


The  UDF  simulation  tool  developed  has  great  flexibility  for  injecting  thermal 
energy  and  momentum  into  the  flow.  Several  variables  maybe  tailored,  as  well  as  the 
weighting  functions  that  control  the  force  density  and  power  density  for  simulation  of  the 
DBD  induced  wall- jet.  Fluent®,  the  commercial  software  used  for  these  simulations, 
was  fast  and  accurate  with  great  capability  and  flexibility.  The  addition  of  the  UDF 
source  terms  to  simulate  a  DBD  was  much  less  complex  and  easier  to  use  than  the 
creation  of  a  special  ionized  flow  solver.  A  plasma  simulation,  constrained  by  time  steps 
of  10’ s  of  nanoseconds  would  take  an  excessive  amount  of  time  and  processing  power. 
The  code  presented  mitigates  this  problem  by  simplifying  the  issue  to  thermal  energy  and 
momentum  addition  into  the  system  using  temporal  averages  of  these  two  sources.  The 
code  has  the  ability  to  take  into  account  a  varying  time  step  that  can  be  smaller  or  larger 
than  the  driving  voltage  wave’s  period.  Many  other  parameters,  such  as  the  mathematical 
description  of  the  spatial  weighting  of  cells  for  simulation  of  the  DBD,  the  number  of 
DBD  elements,  and  the  driving  voltage  frequency  can  also  be  modified.  However, 
extreme  care  is  necessary  when  implementing  the  UDF  source  term  to  have  it  operate 
properly. 

Given  the  two  force  density  profiles  of  Boeuf  and  Pitchford  and  of  Roy  and 
Gaitonde,  the  profile  that  simulates  the  momentum  addition  into  the  flow  the  best  is  the 
Roy  and  Gaitonde  force  density  profile.  The  limited  geometry  size  of  the  Boeuf  and 
Pitchford  force  density  profile  did  not  allow  the  ability  to  overcome  the  local  flow  at  low 
velocities,  which  is  not  what  is  seen  in  experiments.  Further,  when  the  horizontal  extent 
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of  the  Boeuf  and  Pitchford  force  density  profile  was  increased,  minimal  change  was 
observed  in  the  downstream  velocity  profiles.  These  jets  dissipated  quickly  after  the 
force  propelling  them  was  removed.  The  Roy  and  Gaitonde  force  density  profile  appears 
to  address  the  momentum  addition  correctly  by  having  a  significant  vertical  extent  as 
well  as  horizontal.  Also,  the  thermal  addition  effect  for  Roy  and  Gaitonde  is  consistent 
with  estimates  from  section  III  and  simulation  results  depicted  in  Figure  25  for 
temperature  increases. 

Even  though  the  majority  of  power  is  put  into  the  thermal  source,  it  is  my 
conclusion  that  DBD  performance  is  tied  to  the  momentum  source  as  opposed  to  the 
thermal  source.  This  was  demonstrated  in  the  results  of  Figure  57,  when  the  momentum 
source  was  fixed  at  0.18  W/m  input  power  and  the  thermal  source  varied  by  nearly  4 
orders  of  magnitude  from  0.02  W/m  to  199.82  W/m  with  minimal  wall  jet  velocity 
difference  in  the  cases.  Font  [19]  performed  a  further  study  on  this  effect  using  an  air 
chemistry  model  that  included  both  positive  and  negative  ions  and  came  to  the  same 
conclusion. 

The  DBD  performance  does  not  have  fixed  coefficients;  that  is,  the  percentage  of 
power  going  towards  thermal  addition  and  momentum  addition  are  not  fixed  as  the  power 
is  increased.  If  this  were  the  case,  then  the  velocity  of  the  wall  jet  would  increase  as  the 
voltage,  and  subsequent  power,  was  increased  in  experiments.  This  was  not  the  case  as 
reported  by  Newcamp  [18]  in  experiments  and  shown  in  simulations  previously 
discussed. 
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The  airfoil  simulations  concluded  that  Corke’s  data  [7]  could  be  successfully 
simulated  with  a  Roy  and  Gaitonde  force  density  profile.  The  resulting  comparison 
between  the  Fluent®  simulations  and  Corke’s  data  showed  a  similar  minimal 
performance  improvement  trend  up  to  12°  AoA.  This  capability  should  lead  to  faster, 
lower  fidelity,  simulations  of  DBD’s  on  airfoils  and  low-Rec  turbine  blade  research. 

Suggested  follow-on  research  is  to  modify  for  frequency  response  of  the  system, 
vary  placement  of  the  DBD  on  a  NACA  0009  airfoil  to  examine  system  performance,  and 
to  vary  the  number  of  DBD  devices  and  examine  their  effects.  The  time  steps  taken  were 
only  able  to  resolve  10  cycles  of  the  10  kHz  driving  voltage  to  an  average.  To  examine 
the  effectiveness  of  frequency  and  input  waveform  on  the  system  as  Likhanskii  [20]  did, 
the  time  steps  would  need  to  be  lowered  such  that  several  sample  s  are  taken  for  every 
period  to  accurately  resolve  a  temporal  solution.  For  this  follow-on  effort,  limited 
modification  to  the  thermal  energy  source  and  momentum  source  may  be  required.  DBD 
placement  on  the  NACA  0009  airfoil  will  need  to  compare  against  experimental  data. 

This  data  will  need  to  be  generated  and  compared  with  the  follow-on  simulations  for 
DBD  placement.  The  final  parameter  suggested  to  be  modified  under  follow-on  research 
is  the  number  of  DBD’s  and  their  spacing.  The  code  is  set  up  to  add  several  areas  of 
weighting  for  a  multiple  source  DBD  simulation  at  a  fixed  interval  distance.  All  sources 
would  behave  identically  with  respect  to  input  power  to  the  flow,  momentum  addition, 
and  thermal  energy  addition.  Effects  of  increased  lift  while  maintaining  drag  should  be 
examined  and  compared  with  Corke’s  data  [7]  as  the  number  of  simulated  DBD’s 
increases. 
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The  UDF  simulation  tool  developed  has  great  flexibility  for  injecting  thermal 
energy  and  momentum  into  the  flow  with  follow-on  research  that  can  be  done  from  the 
existing  code.  The  creation  of  the  UDF  code  allowed  for  examination  and  macro¬ 
analysis  of  the  Boeuf  and  Pitchford  and  the  Roy  and  Gaitonde  DBD  force  density  profiles 
on  a  flat  plate  geometry.  The  UDF  code  then  utilized  the  Roy  and  Gaitonde  force  density 
profile  to  simulate  a  DBD  on  a  NACA  0009  airfoil  to  compare  results  with  Corke’s  [7] 
experimental  data.  These  efforts  were  the  purpose  of  this  research. 


V-4 


VI.  Appendix  A 


The  tables  that  follow  contain  the  data  points  used  to  create  the  NACA  0009 
airfoil.  These  points  scale  to  a  cord  length  of  0.202  meters  and  represent  only  the  top 
portion  of  the  airfoil.  The  bottom  portion  is  a  mirror  image  as  this  is  a  symmetric  airfoil. 


Table  10.  Computer  Generated  NACA  0009  Airfoil  Data  Point  Set 


X 

X 

Y 

0.0000000 

0.0000017 

0.0000000 

0.0000785 

0.0000873 

0.0005561 

0.0069337 

0.0045703 

0.0000034 

0.0001109 

0.0000890 

0.0005615 

0.0071003 

0.0046189 

0.0570650 

0.0090821 

0.0000051 

0.0001358 

0.0000907 

0.0005668 

0.0072670 

0.0046667 

0.0585800 

0.0090893 

0.0000068 

0.0001567 

0.0000924 

0.0005721 

0.0074336 

0.0047137 

0.0600950 

0.0090924 

0.0000086 

0.0001752 

0.0000942 

0.0005773 

0.0076003 

0.0047601 

0.0616100 

0.0090917 

0.0000103 

0.0001919 

0.0000959 

0.0005825 

0.0077669 

0.0048058 

0.0631250 

0.0090872 

0.0000120 

0.0002072 

0.0000976 

0.0005876 

0.0079336 

0.0048509 

0.0646400 

0.0090791 

0.0000137 

0.0002214 

0.0000993 

0.0005927 

0.0081002 

0.0048954 

0.0661550 

0.0090674 

0.0000154 

0.0002348 

0.0001010 

0.0005977 

0.0082669 

0.0049392 

0.0676700 

0.0090523 

0.0000171 

0.0002475 

0.0001010 

0.0005977 

0.0084335 

0.0049824 

0.0691850 

0.0090338 

0.0000188 

0.0002595 

0.0002677 

0.0009672 

0.0086002 

0.0050251 

0.0707000 

0.0090122 

0.0000205 

0.0002710 

0.0004343 

0.0012266 

0.0087668 

0.0050672 

0.0722150 

0.0089873 

0.0000223 

0.0002820 

0.0006010 

0.0014377 

0.0089335 

0.0051087 

0.0737300 

0.0089595 

0.0000240 

0.0002926 

0.0007676 

0.0016197 

0.0091001 

0.0051497 

0.0752450 

0.0089286 

0.0000257 

0.0003028 

0.0009343 

0.0017817 

0.0092668 

0.0051902 

0.0767600 

0.0088949 

0.0000274 

0.0003127 

0.0011009 

0.0019290 

0.0094334 

0.0052302 

0.0782750 

0.0088584 

0.0000291 

0.0003223 

0.0012676 

0.0020647 

0.0096001 

0.0052697 

0.0797900 

0.0088192 

0.0000308 

0.0003316 

0.0014342 

0.0021911 

0.0097667 

0.0053087 

0.0813050 

0.0087773 

0.0000325 

0.0003406 

0.0016009 

0.0023098 

0.0099334 

0.0053473 

0.0828200 

0.0087328 

0.0000342 

0.0003494 

0.0017675 

0.0024219 

0.0101000 

0.0053853 

0.0843350 

0.0086859 

0.0000359 

0.0003580 

0.0019342 

0.0025283 

0.0116150 

0.0057121 

0.0858500 

0.0086364 

0.0000377 

0.0003664 

0.0021008 

0.0026297 

0.0131300 

0.0060081 

0.0873650 

0.0085846 

0.0000394 

0.0003746 

0.0022675 

0.0027268 

0.0146450 

0.0062780 

0.0888800 

0.0085305 

0.0000411 

0.0003826 

0.0024341 

0.0028200 

0.0161600 

0.0065255 

0.0903950 

0.0084741 

0.0000428 

0.0003904 

0.0026008 

0.0029096 

0.0176750 

0.0067531 

0.0919100 

0.0084155 

0.0000445 

0.0003981 

0.0027674 

0.0029960 

0.0191900 

0.0069632 

0.0934250 

0.0083547 

0.0000462 

0.0004056 

0.0029341 

0.0030796 

0.0207050 

0.0071575 

0.0949400 

0.0082919 

0.0000479 

0.0004130 

0.0031007 

0.0031605 

0.0222200 

0.0073374 

0.0964550 

0.0082270 

0.0000496 

0.0004203 

0.0032674 

0.0032389 

0.0237350 

0.0075042 

0.0979700 

0.0081601 

0.0000514 

0.0004274 

0.0034340 

0.0033150 

0.0252500 

0.0076589 

0.0994850 

0.0080912 

0.0000531 

0.0004344 

0.0036007 

0.0033890 

0.0267650 

0.0078024 

0.1010000 

0.0080204 

0.0000548 

0.0004413 

0.0037673 

0.0034611 

0.0282800 

0.0079355 

0.1058095 

0.0077836 

0.0000565 

0.0004481 

0.0039340 

0.0035313 

0.0297950 

0.0080589 

0.1106190 

0.0075293 

0.0000582 

0.0004548 

0.0041006 

0.0035998 

0.0313100 

0.0081731 

0.1154286 

0.0072586 

0.0000599 

0.0004614 

0.0042673 

0.0036666 

0.0328250 

0.0082786 

0.1202381 

0.0069725 

0.0000616 

0.0004679 

0.0044339 

0.0037319 

0.0343400 

0.0083760 

0.1250476 

0.0066719 

0.0000633 

0.0004743 

0.0046006 

0.0037958 

0.0358550 

0.0084658 

0.1298571 

0.0063576 

0.0000651 

0.0004806 

0.0047672 

0.0038582 

0.0373700 

0.0085482 

0.1346667 

0.0060302 

0.0000668 

0.0004869 

0.0049339 

0.0039194 

0.0388850 

0.0086236 

0.1394762 

0.0056902 

0.0000685 

0.0004930 

0.0051005 

0.0039793 

0.0404000 

0.0086924 

0.1442857 

0.0053381 

0.0000702 

0.0004991 

0.0052672 

0.0040381 

0.0419150 

0.0087548 

0.1490952 

0.0049742 

0.0000719 

0.0005051 

0.0054338 

0.0040957 

0.0434300 

0.0088112 

0.1539048 

0.0045987 

0.0000736 

0.0005110 

0.0056005 

0.0041522 

0.0449450 

0.0088617 

0.1587143 

0.0042117 

0.0000753 

0.0005169 

0.0057671 

0.0042077 

0.0464600 

0.0089067 

0.1635238 

0.0038133 

0.0000770 

0.0005227 

0.0059338 

0.0042622 

0.0479750 

0.0089463 

0.1683333 

0.0034034 

0.0000787 

0.0005284 

0.0061004 

0.0043157 

0.0494900 

0.0089808 

0.1731429 

0.0029820 

0.0000805 

0.0005341 

0.0062671 

0.0043683 

0.0510050 

0.0090103 

0.1779524 

0.0025487 

0.0000822 

0.0005397 

0.0064337 

0.0044201 

0.0525200 

0.0090350 

0.1827619 

0.0021032 

0.0000839 

0.0005452 

0.0066004 

0.0044710 

0.0540350 

0.0090551 

0.1875714 

0.0016452 

0.0000856 

0.0005507 

0.0067670 

0.0045211 

0.0555500 

0.0090707 

0.1923810 

0.2020000 

0.0011742 

0.0000000 
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Table  11.  Data  Set  #1  [12] 


Table  12.  Data  Set  #2  [14] 


X 

y 

1.000000 

0.000000 

0.992588 

0.000977 

0.978431 

0.002727 

0.962932 

0.004519 

0.946984 

0.006252 

0.930890 

0.007916 

0.914755 

0.009526 

0.898598 

0.01 1088 

0.882422 

0.012606 

0.866228 

0.014086 

0.850019 

0.015533 

0.833798 

0.016950 

0.817571 

0.018340 

0.801341 

0.019703 

0.785109 

0.021040 

0.768877 

0.022350 

0.752645 

0.023634 

0.736410 

0.024889 

0.720173 

0.026117 

0.703932 

0.027318 

0.687689 

0.028491 

0.671446 

0.029638 

0.655203 

0.030756 

0.638962 

0.031844 

0.622721 

0.032902 

0.606482 

0.033927 

0.590242 

0.034920 

0.574002 

0.035878 

0.557763 

0.036802 

0.541527 

0.037689 

0.525293 

0.038538 

0.509062 

0.039347 

0.492833 

0.040113 

0.476607 

0.040835 

0.460383 

0.041511 

0.444162 

0.042138 

0.427944 

0.042716 

0.411732 

0.043241 

0.395527 

0.043709 

0.379330 

0.044117 

0.363143 

0.044460 

0.346965 

0.044734 

0.330794 

0.044937 

0.314630 

0.045063 

0.298486 

0.045109 

0.282373 

0.045063 

0.266276 

0.044915 

0.250174 

0.044663 

0.234093 

0.044308 

0.218057 

0.043835 

0.202046 

0.043231 

0.186062 

0.042490 

0.170138 

0.041598 

0.154273 

0.040535 

0.138473 

0.039283 

0.122778 

0.037819 

0.107202 

0.036110 

0.091745 

0.034125 

0.076471 

0.031840 

0.061545 

0.029225 

0.047280 

0.026246 

0.034392 

0.022919 

0.023705 

0.019346 

0.015757 

0.015817 

0.010145 

0.012508 

0.006356 

0.009616 

0.003720 

0.007082 

0.001851 

0.004763 

0.000663 

0.002678 

0.000076 

0.000856 

0.000000 

0.000000 

X 

y 

1.00000 

0.00000 

0.99572 

0.00057 

0.98296 

0.00218 

0.96194 

0.00463 

0.93301 

0.00770 

0.89668 

0.01127 

0.85355 

0.01522 

0.80438 

0.01945 

0.75000 

0.02384 

0.69134 

0.02823 

0.62941 

0.03247 

0.56526 

0.03638 

0.50000 

0.03978 

0.43474 

0.04248 

0.37059 

0.04431 

0.33928 

0.04484 

0.30866 

0.04509 

0.27886 

0.04504 

0.25000 

0.04466 

0.22221 

0.04397 

0.19562 

0.04295 

0.17033 

0.04161 

0.14645 

0.03994 

0.12408 

0.03795 

0.10332 

0.03564 

0.08427 

0.03305 

0.06699 

0.03023 

0.05156 

0.02720 

0.03806 

0.02395 

0.02653 

0.02039 

0.01704 

0.01646 

0.00961 

0.01214 

0.00428 

0.00767 

0.00107 

0.00349 

0.00000 

0.00000 

Figure  67 

Airfoil  Data  Set  #1  Plot  of  Table  1 1 


Figure  68 

Airfoil  Data  Set  #1  Plot  of  Table  12 
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AoA  vs  CL 


Re=1 50600  -g-Re=199900  -o-Re=302900  -^-SATurb  Re=1 80000  -^SATurb  Re=360000 


Figure  69 

NACA  0009  Angle  of  Attack  VS  Lift  Coefficient 
for  Spalart-Allmaras  Turbulence  Model 


AoA  vs  CL 


SATurb  Re=1 80000  — SATurb  Re=360000 


Figure  70 

NACA  0009  Angle  of  Attack  VS  Lift  Coefficient  Difference 
for  Spalart-Allmaras  Turbulence  Model  as  Compared  to  Experimental  Data  by  Selig  [15] 
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AoA  vs  CD 


-o-  Re=1 50600  -g-  Re=199900  -a-  Re=302900  SATurb  Re=180000  -^SATurb  Re=360000 

Figure  7 1 

NACA  0009  Lift  Coefficient  VS  Drag  Coefficient 
for  Spalart-Allmaras  Turbulence  Model 


AoA  vs  CD 


— *-  SATurb  Re=1 80000  SATurb  Re=360000 


Figure  72 

NACA  0009  Lift  Coefficient  VS  Drag  Coefficient  Difference 
for  Spalart-Allmaras  Turbulence  Model  as  Compared  to  Experimental  Data  by  Selig  [15] 
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VII.  Appendix  B 


The  appendix  describes,  in  detail,  the  workings  of  each  part  of  the  Fluent® 
compiled  C  code  and  how  verification  of  each  piece  was  performed.  The  code  has 
several  areas  that  had  to  be  commented  to  allow  compilation  by  Fluent’ s®  C  compiler. 

The  file  that  contains  all  of  the  following  code  is 
“temp_mom_src_trm_FLUENT.c”  and  has  only  one  header  file  for  its  function 
definitions  that  Fluent®  requires.  The  code  is  detailed  in  the  order  it  is  written  in  the  file. 
The  most  recent  version  of  the  code  was  modified  on  19  August  2005. 

Include  Files 

Originally,  the  C++  include  files  were  used  in  this  program.  However  the 

Fluent®  C  compiler  was  unable  to  identify  them,  and  as  a  result,  the  standard  C  files  had 

to  be  included  instead  via  the  udf.h  include  file. 

/*  Fluent  include  files*/ 

#include  "udf.h" 
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Function  Definitions 


Function  definitions  are  required  at  the  beginning  of  every  C  or  C++  program  to 
identify  the  subroutines  that  are  in  the  following  code.  The  asterisks  tell  the  compiler  that 
a  pointer  is  to  be  used.  This  is  necessary  because  several  variables  are  sometimes  needed 
to  be  modified  and  returned  to  the  previous  subroutine. 

Each  subroutine  will  be  described  in  detail  as  to  its  inner  workings  and 
verification  method. 


/* 

Subroutine:  Function  Definitions 


Modified: 
Created : 
Creator : 


04/07/2005 

04/04/2005 

Capt  Timothy  R.  Klein 


Description:  Defines  Functions  to  be  used  throughout  file  in  all 
subroutines . 

*/ 


void  read_data  (); 

real  power_avg  (real, real); 
real  power_f unct ( real ) ; 

void  coord_xform ( real , real ,  real  *,real  *,  real, real,  real, real); 

real  four_point  (real,  real,  real,  real,  real, real,  real, real,  real); 
real  three_point ( real , real ,  real, real,  real, real,  real); 

real  line_side  (  int , real , real ) ; 

void  line_intercept (int,  real,  real,  real,  real,  real  *,real  *); 
real  line_of f set ( int ) ; 


real  curve_y  (real, real); 

real  curve_dy  (real, real); 

real  volume_integration (real,  real,  real, real,  real, real,  real); 
real  weight_funct (real ,  real )  ; 
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Constant  Definitions 


The  research  required  that  certain  parameters  be  easily  changed  from  one  run  of 
data  to  the  next.  Therefore,  it  was  necessary  to  create  a  global  variable  set  as  written 
below. 

/* 

'k'k'k'k'k'k'k'k'k'k'k'k'ki<'k'k'ki<'k'k'ki<'k'k'ki<'k'k'ki('k'k'ki<'k'k'k'k'k'k'ki<'k'k'ki<'k'k'k'k'k'k'k'k'k'k'ki<'k'k'ki<'k'k'k'k'k'k'k'k'k 

Subroutine:  Constant  Definitions 

Modified:  08/17/2005 

Created:  03/21/2005 

Creator:  Capt  Timothy  R.  Klein 

Description:  Defines  Constants  to  be  used  throughout  file  in  all 
subroutines . 

*/ 

/*  Boltzmann  constant,  Joules/Kelvin  */ 

#def ine  kb  1.3807E-23 

/*  Total  Power  in  Watts  per  Unit  Length  (1.0  meter)  */ 
tdefine  POWER_TOT  80.0 

/*  TEMP_PERCENT  +  MOM_PERCENT  =  100.0%  */ 

#def ine  TEMP_PERCENT  98.0 

#def ine  MOM_PERCENT  2.0 

/*  X_MOM  +  Y_MOM  =  100.0%  of  MOM_PERCENT  */ 
tdefine  X_MOM  99.0 

tdefine  Y_MOM  1.0 

/*  Number  of  DBD's*/ 

#def ine  DBD_NUM  1 

/*  Distance  in  meters  between  DBD's  if  multiple  */ 

#def ine  DBD_S PACING  0.050 

/*  Physical  Location  of  first  DBD  (Source) 

If  multiple  DBD's,  then  X_POSITION*CORD  +  n*DBD_SPACING 
Where  n  is  the  number  of  the  DBD  starting  at  0 
Distance  based  on  a  1  meter  cord  length  */ 
tdefine  X_POSITION  0.75  /*  cord  */ 

#def ine  Y_POSITION  0.0  /*  cord  */ 

/*  Frequency  of  DBD  Voltage  in  Hz  */ 
tdefine  FREQ  10000.0 
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/*  Cord  Length  (Flat  Plate  or  Airfoil)  */ 

#define  CORD  0.202  /*  meter  */ 

/*  Airfoil  Span  for  3D  Calculations*/ 

#define  SPAN  1.0  /*  meter  */ 

/*  Switch  for  Flat-Plate  or  Airfoil 
Flat-Plate:  F_OR_A  =  0 
Airfoil:  F_OR_A  =  1  */ 

tdefine  F_OR_A  0 

/*  Switch  for  Bouef f SPitchf ord  or  Roy&Gaitonde  Weighting-Scheme 
Bouef f SPitchford:  WScheme  =  1 
Roy&Gaitonde:  WScheme  =  2  */ 

tdefine  WScheme  2 
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Cell  Weighting  Define  on  Demand  UDF  Code 

This  section  of  code  defines  the  Cell  Weighting  Define  on  Demand  UDF  that  was 
used  in  Fluent®. 


/* 

Subroutine:  DEFINE_ON_DEMAND (cell_weight_on_demand) 


Modified: 
Created : 
Creator : 


08/08/2005 

03/21/2005 

Capt  Timothy  R.  Klein 


Description:  Loops  over  all  cells  to  weight  each  one.  Total  Area  is 
found  and  then  normalized  by  the  Total  Area  to  give  1. 

*/ 

DEF INE_ON_DEMAND ( ce 1 l_we ight_on_demand) 


/*  Define  Fluent  cell  variables 
Domain  *d; 

Thread  *t; 
cell_t  c; 

Node  *node; 


d  =  Get_Domain  ( 1 ) ; 
t  =  Lookup_Thread (d,  2); 


/*  Get  fluid  thread  using  Fluent  utility  */ 
/*  Zone=2  is  for  the  fluid  */ 


real  x  [4]  ,  y  [4]  ; 
real  u [ 4 ] ,  v  [  4 ]  ; 

int  tot_nodes,  count; 
int  i,  n; 


real  tot_weight  =0.0;  / i 

real  weight  =  0.0;  / i 


Set  initial  Total  Weight  */ 

Set  initial  Weight  for  the  cell  */ 


if  (N_UDM< (DBD_NUM+3+l) ) 


printf ("\n\n\nYOU  MUST  DEFINE  %d  UDFM ' s !  ! \n"  ,  (DBD_NUM+3  +  l )  )  ; 
Internal_Error ( "YOU  MUST  DEFINE  more  UDFM ' s !  ! \n\n\n " )  ; 
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/*  Fill  the  UDM  (User  Defined  Memory)  with  cell  weight  */ 
thread_loop_c (t,  d) 

{ 


/*  Loop  for  Multiple  DBD ' s  */ 
for (i=l ; i<=DBD_NUM; ++i) 

{ 

/*  Set  Location  of  DBD  and  find  Parallel  Vector  for  Coordinate 
Transformations  */ 

real  px  =  CORD*X_POSITION  +  ( (real) ( 1—1 ) *DBD_SPACING) ; 

real  py  =  curve_y (px, CORD) ; 

real  vx  =  1.0; 

real  vy  =  curve_dy (px, CORD ) ; 

real  vmag  =  sqrt (vx*vx+vy*vy ) ; 

vx  /=  vmag; 
vy  /=  vmag; 

/*  Weight  each  cell  */ 
begin_c_loop (c, t )  /*  Cell  Loop  */ 

{ 

C_UDMI (c, t, i)  =  0.0; 

/*  Get  the  total  number  of  nodes  in  a  cell  */ 
tot_nodes  =  C_NNODES (c, t ) ; 

/*  Get  the  coordinates  of  the  nodes  in  a  cell  */ 
c_node_loop (c, t, count)  /*  Node  Loop  */ 

{ 

node  =  C_NODE (c, t, count) ; 
x[count]  =  NODE_X (node ) ; 
y[ count]  =  NODE_Y (node) ; 

} 

count  =  0; 

weight  =0.0;  /*  Set  initial  Weight  for  the  cell  */ 

for ( count =0 ; count <t ot_node s ; t+count ) 

{ 

coord_xform(x [count] ,y [count] ,  &u [count ], &v [count ] ,  vx,vy. 
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/*  Send  points  to  functions  to  see  if  a  cell  straddles  a  line 
and  to  give  weight  */ 
if (tot_nodes  ==  4) 

{ 

weight  =  four_point  (u[0],v[0],  u[l],v[l],  u[2],v[2], 

u[3],v[3],  weight); 

} 

else  if (tot_nodes  ==  3) 

{ 

weight  =  three_point (u [ 0 ] , v [ 0 ] ,  u[l],v[l],  u[2],v[2], 

weight ) ; 

} 

else 

{ 

printf ( "ERROR [ cell_weight_on_demand] :  Total  Number  of  Nodes 
NOT  3  or  4  !  \ n " ) ; 

} 

C_UDMI (c, t , i )  =  weight; 
tot_weight  +=  fabs (weight ) ; 

if  (c%  (100*100) ==0)  printf (" \n" ); /*New  line  every 
if (c%100  ==0)  printf  (".");  /*Shows  progress 

} 

end_c_loop (c, t )  /*  Cell  Loop  */ 

/*  Normalize  each  cell  */ 
begin_c_loop (c, t )  /*  Cell  Loop  */ 

{ 

C_UDMI (c, t, i)  /=  tot_weight;  /*Normalize  each  cell*/ 

if (c% (100*100) ==0)  print f (" \n" );/ *New  line  every  100  "."*/ 

if (c%100  ==0)  printf (".");  /*Shows  progress  of  cell  loop*/ 

} 

end_c_loop (c, t )  /*  Cell  Loop  */ 

}  /*  DBD_NUM  for  loop  */ 

}  /*  Thread  Loop  */ 

}  /*  DEFINE_ON_DEMAND (cell_weight_on_demand)  */ 


10  0  "  .  " * / 
of  cell  loop*/ 
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Temperature  Source  UDF  Code 

This  section  of  code  defines  the  Temperature  Source  UDF  that  was  used  in 
Fluent®. 


/* 

■k-k-kit-kit-kit-k-k-kit-kit-kit-kit-kiic-kit-kit-kit-kit-kit-kit-k-k-kiic-k-k-kiic-k-k-kit-k-k-kit-kit-kit-kit-kit-k-k-kiic-kit-kit-kit-kit-kit-k 

Subroutine :  DEFINE_SOURCE ( temp_source ) 


Modified: 
Created : 
Creator : 


08/19/2005 

03/21/2005 

Capt  Timothy  R.  Klein 


Description:  Adds  thermal  energy  to  cell  given  the  cell,  power,  and 
weighting . 

*/ 

DEF INE_SOURCE (temp_source, c, t, dS, eqn) 

{ 


if (TEMP_PERCENT>0 . 0)  /*  Less  Computational  Time  */ 

{ 

real  weight  =  0.0; 
real  source; 

int  i; 

/*  Get  normalized  cell  weight  from  memory  */ 

/*  i  represents  each  DBD  if  multiple  */ 

for (i=0 ; i<DBD_NUM; ++i) 

{ 

weight  +=  C_UDMI (c, t, i+1) ; 

} 


i f (weight  ! =  0.0) 

{ 

real  volume  =  C_VOLUME (c, t ) ; 


/*  This  needs  to  be  a  volume  ratio  of  W/mA3  */ 

source  =  POWER_TOT* (TEMP_PERCENT/100 . 0) *fabs (weight) /volume; 

C_UDMI ( c, t , DBD_NUM+3 )  =  source;  /*if  DBD_NUM=1  then  UDMI  slot=4*/ 

} 

else 


source  =  0.0; 

C_UDMI (c, t, DBD_NUM+3)  =  source;  /*if  DBD_NUM=1  then  UDMI  slot=4*/ 
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dS [eqn]  =  0.0; 
return  source; 


}  /*  TEMP_PERCENT  >  0.0  */ 

else  /*  TEMP_PERCENT  <=  0.0  */ 

{ 

C_UDMI ( c, t , DBD_NUM+3 )  =  0.0;  /*if  DBD_NUM=1  then  UDMI  slot=4*/ 
dS [eqn]  =  0.0; 
return  0.0; 


}  /*  DEFINE_SOURCE (temp_source)  subroutine  end  */ 
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X-Momentum  Source  UDF  Code 

This  section  of  code  defines  the  X-Momentum  Source  UDF  that  was  used  in 


Fluent®. 


/* 

Subroutine:  DEFINE_SOURCE ( x_raoment  um_s  ource ) 


Modified: 
Created : 
Creator : 


08/19/2005 

03/21/2005 

Capt  Timothy  R.  Klein 


Description:  Adds  x-momentum  energy  to  cell  given  the  cell,  power,  and 
weighting . 

*/ 

DEFINE_SOURCE (x_momentum_source, c, t, dS, eqn) 


if (MOM_PERCENT>0 . 0)  /*  Less  Computational  Time  */ 


real  weight  =  0.0; 
real  source; 
real  u02,  v02; 


int  i; 

/*  Get  normalized  cell  weight  from  memory  */ 
/*  i  represents  each  DBD  if  multiple  */ 
for (i=0 ; i<DBD_NUM; ++i) 


weight  +=  C_UDMI (c, t , i+1) ; 

} 


real  dt  =  CURRENT_TIMESTEP ;  /*  Get  Current  Timestep  */ 


i f (weight  ! =  0.0) 

{ 

real  energy  =  POWER_TOT*dt ; 

/*  Calculate  weighted  energy  into  cell  */ 

real  energy_x  =  (MOM_PERCENT/100 . 0) * (X_MOM/100 . 0) *weight*energy; 
real  energy_y  =  (MOM_PERCENT/100 . 0) * (Y_MOM/100 . 0) *weight*energy; 

/*  Calculate  mass  in  cell  */ 
real  density  =  C_R(c,t); 
real  volume  =  C_VOLUME (c,  t ) ; 
real  mass  =  volume*density; 

/*  Get  Initial  Cell  Velocities  and  Kinetic  Energy  */ 
real  uOO  =  C_U(c,t);  /*  initial  u  velocity  */ 

real  vOO  =  C_V(c,t);  /*  initial  v  velocity  */ 
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/  *  kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk J 

/*  Allow  for  direction  of  source  to  be  along  surface  tangent  */ 


/*  Set  Surface  Tangent  Vector  at  DBD  Location  */ 

real  px  =  CORD*X_POSITION  +  (  (real)  (i) *DBD_SPACING) ; 

real  vx  =  1.0; 

real  vy  =  curve_dy (px, CORD )  ; 


/*  Magnitude  of  vector  v  */ 

real  v_xy_mag  =  sqrt (vx*vx  +  vy*vy) ; 
vx  /=  v_xy_mag; 
vy  /=  v_xy_mag; 

/*Reorient  x&y  velocity  vectors  to  be  tangent  to  surface  at  DBD*/ 
real  uOl  =  u00*vx  +  v00*vy; 
real  vOl  =  -u00*vy  +  v00*vx; 

/*  Calculate  Cell  Velocities  and  Kinetic  Energy  after  energy 
addition  */ 
real  KElx  =  energy_x; 
real  KEly  =  energy_y; 

/*  Need  to  preserve  Energy  direction  by  keeping  velocity  sign  */ 
/*  Avoids  sqrt  of  negative  number  and  allows  for  deceleration  */ 
/*  if  flow  moving  in  opposite  direction  of  DBD  wall  jet  */ 

real  u02sqr  =  u01*fabs (uOl)  +  (2* (KElx) /mass ) ; 
if (u02sqr  <  0.0) 

{ 

u02  =  -  sqrt (  fabs (  u02sqr  )  ) ; 

} 

else 

{ 

u02  =  sqrt  (  u02sqr  )  ; 

} 


/*  Need  to  preserve  Energy  direction  by  keeping  velocity  sign  */ 
/*  Avoids  sqrt  of  negative  number  and  allows  for  deceleration  */ 
/*  if  flow  moving  in  opposite  direction  of  DBD  wall  jet  */ 

real  v02sqr  =  v01*fabs (vOl)  +  (2*  (KEly) /mass) ; 
if (v02sqr  <  0.0) 

{ 

v02  =  -  sqrt (  fabs (  v02sqr  )  ) ; 

} 

else 

{ 

v02  =  sqrt  (  v02sqr  )  ; 

} 

/*  Return  to  x&y  coordinate  system  */ 
real  ull  =  u02*vx  -  v02*vy; 
real  acl  =  (ull-uOO) /dt; 

/ k  kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk  k / 
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/*  Density  instead  of  mass  so  there  is  a  volume  ratio  */ 

/*  momentum  source  =  N/mA3  =  density  *  acceleration  */ 

source  =  density*acl;  /*  Debug  */ 

C_UDMI ( c, t , DBD_NUM+ 1 )  =  source;  /*if  DBD_NUM=1  then  UDMI  slot=2*/ 

} 

else 

{ 

source  =  0.0; 

C_UDMI ( c, t , DBD_NUM+ 1 )  =  source;  /*if  DBD_NUM=1  then  UDMI  slot=2*/ 


dS [eqn]  =  0.0; 
return  source; 


}  /*  MOM_PERCENT  >  0.0  */ 

else  /*  MOM_PERCENT  <=  0.0  */ 

{ 

C_UDMI ( c , t , DBD_NUM+ 1 )  =  0.0;  /*if  DBD_NUM=1  then  UDMI  slot=2*/ 
dS [eqn]  =  0.0; 
return  0.0; 


}  /*  DEFINE_SOURCE (x_momentum_source)  subroutine  end  */ 


VII- 12 


Y-Momentum  Source  UDF  Code 


This  section  of  code  defines  the  Y-Momentum  Source  UDF  that  was  used  in 


Fluent®. 


/* 


Subroutine:  DEFINE_SOURCE ( y_momentum_source ) 


Modified:  08/16/2005 

Created:  03/21/2005 

Creator:  Capt  Timothy  R.  Klein 


Description:  Adds  y-momentum  energy  to  cell  given  the  cell,  power,  and 


weighting . 


*/ 

DEFINE_SOURCE (y_momentum_source, c, t, dS, eqn) 

{ 

if (MOM_PERCENT>0 . 0)  /*  Less  Computational  Time  */ 

{ 

real  source,  u02,  v02,  weight  =  0.0; 
int  i ; 

/*  Get  normalized  cell  weight  from  memory  */ 

/*  i  represents  each  DBD  if  multiple  */ 

for (i  =  0 ;  i<DBD_NUM; +  +  i) 

{ 

weight  =  C_UDMI (c, t , i+1 ) ; 

} 

real  dt  =  CURRENT_TIMESTEP ;  /*  Get  Current  Timestep  */ 

if (weight  ! =  0.0) 

{ 

real  energy  =  POWER_TOT*dt ; 

/*  Keep  the  wall-jet  pulled  to  the  surface  */ 
weight  =  -  fabs (weight ) ; 

/*  Calculate  weighted  energy  into  cell  */ 

real  energy_x  =  (MOM_PERCENT/100 . 0) * (X_MOM/100 . 0) *weight*energy; 
real  energy_y  =  (MOM_PERCENT/100 . 0) * (Y_MOM/100 . 0) *weight*energy; 

/*  Calculate  mass  in  cell  */ 
real  density  =  C_R(c,t); 
real  volume  =  C_VOLUME (c, t ) ; 
real  mass  =  volume*density; 


/*  Get  Initial  Cell  Velocities  and  Kinetic  Energy  */ 
real  uOO  =  C_U(c,t);  /*  initial  u  velocity  */ 
real  vOO  =  C_V(c,t);  /*  initial  v  velocity  */ 


*/ 

*/ 
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/  *  kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk J 

/*  Allow  for  direction  of  source  to  be  along  surface  tangent  */ 


/*  Set  Surface  Tangent  Vector  at  DBD  Location  */ 

real  px  =  CORD*X_POSITION  +  (  (real)  (i) *DBD_SPACING) ; 

real  vx  =  1.0; 

real  vy  =  curve_dy (px, CORD )  ; 


/*  Magnitude  of  vector  v  */ 

real  v_xy_mag  =  sqrt (vx*vx  +  vy*vy) ; 
vx  /=  v_xy_mag; 
vy  /=  v_xy_mag; 

/*Reorient  x&y  velocity  vectors  to  be  tangent  to  surface  at  DBD*/ 
real  uOl  =  u00*vx  +  v00*vy; 
real  vOl  =  -u00*vy  +  v00*vx; 

/*  Calculate  Cell  Velocities  and  Kinetic  Energy  after  energy 
addition  */ 
real  KElx  =  energy_x; 
real  KEly  =  energy_y; 

/*  Need  to  preserve  Energy  direction  by  keeping  velocity  sign  */ 
/*  Avoids  sqrt  of  negative  number  and  allows  for  deceleration  */ 
/*  if  flow  moving  in  opposite  direction  of  DBD  wall  jet  */ 

real  u02sqr  =  u01*fabs (uOl)  +  (2* (KElx) /mass ) ; 
if (u02sqr  <  0.0) 

{ 

u02  =  -  sqrt (  fabs (  u02sqr  )  ) ; 

} 

else 

{ 

u02  =  sqrt  (  u02sqr  )  ; 

} 


/*  Need  to  preserve  Energy  direction  by  keeping  velocity  sign  */ 
/*  Avoids  sqrt  of  negative  number  and  allows  for  deceleration  */ 
/*  if  flow  moving  in  opposite  direction  of  DBD  wall  jet  */ 

real  v02sqr  =  v01*fabs (vOl)  +  (2*  (KEly) /mass) ; 
if (v02sqr  <  0.0) 

{ 

v02  =  -  sqrt (  fabs (  v02sqr  )  ) ; 

} 

else 

{ 

v02  =  sqrt  (  v02sqr  )  ; 

} 

/*  Return  to  x&y  coordinate  system  */ 
real  vll  =  u02*vy  +  v02*vx; 
real  acl  =  (vll-vOO) /dt; 

/ k  kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk/ 
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/*  Density  instead  of  mass  so  there  is  a  volume  ratio  */ 

/*  momentum  =  mass  *  acceleration  */ 

/*  This  requires  momentum_per_volume  =  density  *  acceleration  */ 
/*  Density  is:  mass/volume  */ 

source  =  density*acl;  /*  Debug  */ 

C_UDMI ( c, t , DBD_NUM+2 )  =  source;  /*if  DBD_NUM=1  then  UDMI  slot=3*/ 

} 

else 

{ 

source  =  0.0; 

C_UDMI ( c, t , DBD_NUM+2 )  =  source;  /*if  DBD_NUM=1  then  UDMI  slot=3*/ 

} 


dS [eqn]  =  0.0; 
return  source; 


}  /*  MOM_PERCENT  >  0.0*/ 

else  /*  MOM_PERCENT  <=  0.0  */ 

{ 

C_UDMI ( c, t , DBD_NUM+2 )  =  0.0;  /*if  DBD_NUM=1  then  UDMI  slot=3*/ 
dS [eqn]  =  0.0; 
return  0.0; 


}  /*  DEFINE_SOURCE (y_momentum_source)  subroutine  end  */ 
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Main  Testing  Program  for  Verification 


After  all  of  the  parts  to  this  code  were  finished,  it  was  necessary  to  check  for  their 


correctness  of  content.  The  main  program  that  is  listed  below  does  this  testing,  and  is 
commented  out  for  the  period  of  time  when  the  code  is  inserted  into  Fluent®. 


/* 

'k'k'k'k'ki<'k'k'k'k'k-k'ki<'k'k'k'k'k'k'ki<'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'ki<'k'k'ki<'k'k'ki<'k'k'ki<'k'k'k'k'k'k'k'k'k'k'k'k'k 

Subroutine:  main () 


Modified: 
Created : 
Creator : 


04/08/2005 

04/04/2005 

Capt  Timothy  R.  Klein 


Description:  Tests  different  areas  of  temp_mom_src_trm . cpp  for  compile- 
ability  and  correctness  of  function/subroutine  algorithms. 


*/ 

int  main(int  argc,  char*  argv [ ] ) 


Testing  of  the  settings  is  an  easy  way  to  make  sure  that  you  are  getting  all  of  the  correct 
variables  input. 

//  Test  Settings 

std::cout  <<  "Test  Settings\n\n" ; 
std : : cout 


<< 

"POWER_TOT 

"  << 

It 

20 . 0 

II 

<< 

POWER_TOT 

<< 

"\n 

<< 

" TEMP_PERCENT 

"  << 

II 

100 . 0 

II 

<< 

TEMP_PERCENT 

<< 

"  \  n 

<< 

"MOM_PERCENT 

"  << 

II 

0 . 0 

II 

<< 

MOM_PERCENT 

<< 

"  \  n 

<< 

" X_MOM 

"  << 

II 

90 . 0 

II 

<< 

X_MOM 

<< 

"  \  n 

<< 

" Y_MOM 

"  << 

II 

10 . 0 

II 

<< 

Y_MOM 

<< 

"\n 

<< 

" DBD_NUM 

"  << 

II 

1 

II 

<< 

DBD_NUM 

<< 

"  \  n 

<< 

"DBD_S PACING 

"  << 

II 

0 . 050 

II 

<< 

DBD_S PACING 

<< 

"  \  n 

<< 

" X_POS I T ION 

"  << 

II 

0 . 75 

II 

<< 

X_POSITION 

<< 

"\n 

<< 

" Y_POS I T ION 

"  << 

II 

0 . 0 

II 

<< 

Y_POSITION 

<< 

"  \  n 

<< 

"FREQ 

"  << 

II 

10000 . 0 

II 

<< 

FREQ 

<< 

"  \  n 

<< 

"CORD 

"  << 

II 

0 .202 

II 

<< 

CORD 

<< 

"\n 

<< 

"SPAN 

"  << 

II 

1 . 0 

II 

<< 

SPAN 

<< 

"  \  n 

<<  " F_OR_A 
: cout  <<  "\n\n\n1 

"  << 

. 

I 

II 

0 

II 

<< 

F_OR_A 

<< 

"  \  n 
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The  power  average  and  power  function  subroutines  integrate  the  area  under  the  curve  of 
the  power  function.  Thus,  it  was  tested  that  between  t=0  and  t=7t  radians  that  the  function 
would  return  the  same  value  as  between  t=27t  and  t=37t .  Other  areas  were  compared  to 
ensure  accuracy  of  the  integration  under  the  curve. 

power_avg 


//  Test  subroutine: 
std::cout  <<  "Test  subroutine: 

0.0; 

0 .5/FREQ; 

1 . 0/FREQ; 

1 . 5/FREQ; 

2 . 0/FREQ; 

0 . 3/FREQ; 

1 .3/FREQ; 


&  power_funct 

power_avg  &  power_funct\n\n" ; 


real  tl 
real  t2 
real  t3 
real 
real 
real 
real 


1 4 
t5 
1 6 
1 7 


std : 

: cout  << 

power_avg (t 1 , t2 ) 

<< 

II  II 

<< 

power_avg (t2 , t3 )  << 

"\n" 

<< 

power_avg ( 1 3 ,  t4) 

<< 

II  II 

<< 

power_avg ( 1 4  ,  t5 )  << 

"\n" 

<< 

power_avg (tl,  t3) 

<< 

II  II 

<< 

power_avg (t2, t4 )  << 

II  II 

<< 

power_avg (t3,  t5) 

<< 

"\n" 

// 

should  be  the  same 

as  power. 

_avg ( [ 1 1 , 1 3 ] ,  [ 1 2 , 1 4 ] , 

[t3, t5 

<< 

power_avg ( 1 6 ,  t7) 

<< 

"\n" 

<< 

power_avg (tl,  t5) 

<< 

"\n" ; 

std : 

: cout  << 

"\n\n\n" ; 

Hand  calculations  were  performed  for  a  set  of  points  in  free  space  being  rotated  to  a  new 
coordinate  system  by  a  given  vector.  These  calculations  were  coded  below  and  the 
subroutine  was  sent  the  data  to  give  its  results.  The  program  generated  the  data  expected. 

//  Test  subroutine:  coord_xform 

std::cout  <<  "Test  subroutine:  coord_xf orm\n\n" ; 
real  u,v; 

coord_xf orm ( +  3 . 0 , -2 . 0 ,  u,v,  +2.0, +1.0,  +2.0, +2.0); 

std : : cout  <<  " (2 . 68328, -3 . 1305) \n" ; 
std::cout  <<  "("  <<  u  <<  <<  v  <<  ")\n"; 

std::cout  <<  "\n\n\n"; 


Hand  calculations  were  performed  for  a  set  of  points  in  free  space  on  either  side  of  the 
following  lines.  These  calculations  were  coded  below  and  the  subroutine  was  sent  the 
data  to  give  its  results.  The  program  generated  the  data  expected. 


// 

Test 

subroutine : 

line_of f set 

std 

: : cout 

<< 

"Test  subroutine:  line. 

_of f set\n\n" ; 

std 

: : cout 

<< 

"Last  Column  should  be 

0. 

0  0  0  \  n  " 

<< 

II 

0 

"  << 

line_of f set ( 1 ) 

<< 

0 . 000-line_ 

.offset 

(1) 

<< 

"  \  n  " 

<< 

II 

-0  . 

001  "  << 

line_of f set (2 ) 

<< 

-0 . 001-line_ 

.offset 

(2) 

<< 

"  \  n  " 

<< 

II 

0  . 

01  "  << 

line_of f set ( 3 ) 

<< 

0 . 010-line_ 

.offset 

(3) 

<< 

"  \  n  " 

<< 

II 

0 

"  << 

line_of f set (4 ) 

<< 

0 . 000-line_ 

.offset 

(4) 

<< 

"  \  n  " 

<< 

II 

0  . 

01  "  << 

line_of f set ( 5 ) 

<< 

0 . 010-line_ 

.offset 

(5) 

<< 

"  \  n  "  ; 

std 

: : cout 

<< 

"  \  n  \  n  \  n 

II 
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Hand  calculations  were  performed  for  the  point  (0.005,0.005).  The  subroutine  line_side 
returns  a  double  or  real  value  of  the  distance  from  the  line  number  the  point  is  compared 
against.  These  calculations  were  coded  below  and  the  subroutine  was  sent  the  data  to 
give  its  results.  The  program  generated  the  data  expected. 

//  Test  subroutine:  line_side 

std::cout  <<  "Test  subroutine:  line_side\n\n" ; 
std::cout  <<  line_side ( 1 , +0 . 005 , +0 . 005 )  <<  "  " 

<<  line_side (2, +0 . 005, +0 . 005)  <<  "  " 

<<  line_side (3, +0 . 005, +0 . 005)  <<  "  " 

<<  line_side (4, +0 . 005, +0 . 005)  <<  "  " 

<<  line_side (5, +0 . 005, +0 . 005)  <<  "\n"; 
std::cout  <<  "\n\n\n"; 

Hand  calculations  were  performed  for  the  line  created  by  connecting  the  points 
(0.002,0.003)  and  (-0.002,0.001)  across  line  number  1,  where  x=0.0.  The  returned  values 
of  x  and  y  are  the  location  at  which  the  point-  line  intersects  line  number  1 .  These 
calculations  were  coded  below  and  the  subroutine  was  sent  the  data  to  give  its  results. 

The  program  generated  the  exact  data  expected. 

The  case  where  the  points  do  not  sit  on  either  side  of  the  line  was  also  performed.  An 
error  message  was  built  in  and  successfully  tested  if  this  condition  should  occur. 

//  Test  subroutine:  line_intercept 

std::cout  <<  "Test  subroutine:  line_intercept\n\n" ; 
real  x,y; 

line_intercept ( 1 ,  +0 . 002, +0 . 003,  -0 . 002, +0 . 001,  x,y); 
std::cout  <<  "  (0, 0 . 002) \n" ; 

std::cout  <<  "("  <<  x  <<  <<  y  <<  ")\n\n"; 

line_intercept (2,  +0.002,10.003,  +0.002,10.001,  x,y); 
std::cout  <<  "Point  does  not  stradle  line,  therefore  (0,0) \n"; 
std::cout  <<  "("  <<  x  <<  <<  y  <<  ")\n\n"; 

std::cout  <<  "\n\n\n"; 
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Hand  calculations  were  performed  for  several  points.  The  subroutine  weight_funct 
returns  a  double  or  real  value  from  the  equation  specified  at  the  point  given.  These 
calculations  were  coded  below  and  the  subroutine  was  sent  the  data  to  give  its  results. 
The  program  generated  the  data  expected. 


//  Test  sub 
std : : cout  < 
real  weight 
weight  =  we 
std : : cout  < 
< 

weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
weight  =  we 
std : : cout  < 
< 

weight  =  we 
std : : cout  < 
< 

weight  =  we 
std : : cout  < 
< 

weight  =  we 
std : : cout  < 
< 

std : : cout  < 


routine:  weight_funct 
<  "Test  subrouti 


ight_ 

<  "E : 

<  "7 
ight_ 

<  "0 
ight_ 
<  "0 
ight_ 
<  "0 
ight_ 
<  "0 
ight_ 
<  "0 
ight_ 
<  "0 
ight_ 
<  "0 
ight_ 
<  "0 
ight_ 

<  "E : 

<  "0  . 
ight_ 

<  "E : 

<  "0  . 
ight_ 

<  "E  : 

<  "6. 
ight_ 

<  "E : 

<  "6. 
<  "\n 


0\n" 


-0 . 015\n" 
>851193960: 
:unct  ( —  0 . 0  ( 
-0 . 015\n" 
>851193960: 
:unct  (  +  0 . 0( 
-0 . 0  0  2 \ n " 
>860139906' 
:unct (+0 . 0( 
-0 . 0  0  2 \ n " 
>860139906' 
\n\n" ; 


weight_f unct 

,\n\n" ; 

0 

. 000) ; 

II 

\n" ; 

-0 

.001) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

+  0 

.001) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

+  0 

.011) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

-0 

.001) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

+  0 

.001) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

+  0 

•Oil) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

+  0 

•Oil) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

+  0 

.011) ; 

// 

Outside 

boundary 

area 

II 

\n" ; 

-0 

.001) ; 

// 

Outside 

boundary 

area 

<< 

weight 

<< 

"  \  n  "  ; 

+  0 

.001) ; 

<< 

weight 

<< 

"\n" ; 

-0 

.001) ; 

<< 

weight 

<< 

"\n" ; 

+  0 

.001) ; 

<< 

weight 

<< 

"\n" ; 

Hand  calculations  were  performed  for  a  set  of  3  points.  The  subroutine  returned  the 
weight  of  the  cell  obtained  via  descritized  integration.  These  calculations  were  coded 
below  and  the  subroutine  was  sent  the  data  to  give  its  results.  The  program  generated  the 
data  to  within  0.1%  of  the  expected  value.  Because  this  type  of  integration  is  not  exact 
and  the  volume  step  size  was  set  low  for  speed,  this  is  a  reasonable  answer. 

//  Test  subroutine:  volume_integration 

std:: cout  <<  "Test  subroutine:  volume_integration\n\n" ; 
real  weight; 
weight  = 

volume_integrat ion ( 0 . 00 1 , 0 . 004  ,  0.005,0.008,  0.007,0.001,  0.0); 
std::cout  <<  "0.000125539  "  <<  weight  <<  "\n"; 

std::cout  <<  "\n\n\n"; 
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Hand  calculations  were  performed  for  several  points.  The  subroutine  returned  the  y- 
coordinate  of  the  airfoil  curve,  given  the  airfoil  is  symmetric.  These  calculations  were 
coded  below  and  the  subroutine  was  sent  the  data  to  give  its  results.  The  program 
generated  the  data  to  within  10"6%  of  the  expected  value;  well  within  acceptable  limits. 

//  Test  subroutine:  curve_y 

std::cout  <<  "Test  subroutine:  curve_y\n\n" ; 

real  x; 


x  =  0.000; 
std : : cout  << 
<< 

x  =  0.0101; 
std : : cout  << 
<< 

x  =  0.0707; 
std : : cout  << 
<< 

x  =  0.101; 
std : : cout  << 
<< 

x  =  0.202; 
std : : cout  << 
<< 

std : : cout  << 


" (0, 0) \n" 


("  <<  x  <<  <<  curve_y (x, CORD ) 

<< 

") \n\n"; 

(0 . 0101, 0 . 00538530) \n" 

("  <<  x  <<  <<  curve_y (x, CORD ) 

<< 

") \n\n"; 

(0 . 0707, 0 . 00901220) \n" 

("  <<  x  <<  <<  curve_y (x, CORD ) 

<< 

") \n\n"; 

(0.101,0. 00802040) \n" 

("  <<  x  <<  <<  curve_y (x, CORD ) 

<< 

") \n\n"; 

(0.202,-0) \n" 

("  <<  x  <<  <<  curve_y (x, CORD ) 

<< 

") \n"; 

"\n\n\n" ; 


Hand  calculations  were  performed  for  several  points.  The  subroutine  returned  the  slope 
of  the  airfoil  curve  at  the  x- coordinate,  given  the  airfoil  is  symmetric.  These  calculations 
were  coded  below  and  the  subroutine  was  sent  the  data  to  give  its  results.  The  program 
generated  the  data  to  within  0.1%  of  the  expected  value;  well  within  acceptable  limits. 

//  Test  subroutine:  curve_dy 

std::cout  <<  "Test  subroutine:  curve_dy\n\n" ; 
real  x, slope_data; 
x  =  0.0000000000000001; 

std:: cout  <<  "As  x->0.0,  then  dy->inf inity\n"  ; 
std::cout  <<  " (le-016, infinity) \n" 

<<  "("  <<  x  <<  <<  curve_dy (x, CORD)  <<  ")\n\n"; 

x  =  0.060095; 

slope_data  =  (l./3.)*(  ( 0 . 00 908 93-0 . 00 90 924 )/( 0 . 058580-0 . 0 600 95 )  + 

(0.0090 92 4 -0.0090917)/(0. 060095-0. 061610)  + 

(0. 0090893-0. 0090917)/(0. 058580-0. 061610)  ); 

std::cout  <<  "As  y->max,  then  dy->0.0\n"; 
std : : cout 

<<  "  (0.060095,-0) \n" 

<<  "(0.060095,"  <<  slope_data  <<  ")\n" 


<< 

"  (0  . 

.060095, " 

<< 

(0. 

.  0090893 

-0 

.  0090924) 

/ 

(0. 

,  058580- 

-0  . 

.  0 

60095) 

<< 

" )  \n" 

<< 

"  (0  . 

.060095, " 

<< 

(0. 

.  0090924 

-0 

.  0090917) 

/ 

(0. 

,060095- 

-0  . 

.  0 

61610) 

<< 

" )  \n" 

<< 

"  (0  . 

.060095, " 

<< 

(0. 

.  0090893 

-0 

.  0090917) 

/ 

(0. 

,  058580- 

-0  . 

.  0 

61610) 

<< 

" )  \n" 

<< 

IT  ^  II 

<<  x  << 

II  II 

r 

<< 

curve_dy ( 

x , CORD )  < 

< 

") 

\n\n" ; 

x  =  0.0707; 

slope_data  =  (l./3.)*(  (0 . 0090338-0 . 0090122) / (0 . 069185-0 . 070700)  + 

(0. 0090122-0. 0089873)/(0. 070700-0. 072215)  + 
(0.0090338-0.0089873)/ (0.069185-0 .072215)  )  ; 

std:: cout  <<  "Use  Data  around  x=0.0707  to  get  slope  and  compare\n"; 
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std::cout  <<  "(0.0707,"  <<  slope_data  <<  ")\n" 

<<  "("  <<  x  <<  <<  curve_dy (x, CORD)  <<  ")\n\n"; 

std::cout  <<  "\n\n\n"; 

Hand  calculations  were  performed  for  two  sets  of  3  points.  The  subroutine  returned  the 
weight  of  the  area  contained  within  these  3  points.  These  calculations  were  coded  below 
and  the  subroutine  was  sent  the  data  to  give  its  results.  The  program  generated  the  data  to 
within  0.1%  of  the  expected  value;  well  within  acceptable  limits. 

//  Test  subroutine:  three_point 

std::cout  <<  "Test  subroutine:  three_point\n\n" ; 
real  weight; 

std::cout  <<  "Test  #1:  All  3  points  in  Quadrant  I\n"; 

weight  =  three_point (  0.001,0.004,  0.005,0.008,  0.007,0.001,  0.0); 

std::cout  <<  "0.000125539  "  <<  weight  <<  "\n\n"; 

std : : cout 

<<  "Test  #2:  2  points  in  Quadrant  I,  1  point  in  Quadrant  III\n"; 
weight  =  three_point ( -0 . 00 1 , 0 . 004  ,  0.005,0.008,  0.007,0.001,  0.0); 
std : : cout  <<  "0.000172459176  "  <<  weight  <<  "\n\n"; 

std::cout  <<  "\n\n\n"; 

Hand  calculations  were  performed  for  two  sets  of  4  points.  The  subroutine  returned  the 
weight  of  the  area  contained  within  these  4  points.  These  calculations  were  coded  below 
and  the  subroutine  was  sent  the  data  to  give  its  results.  The  program  generated  the  data  to 
within  0.1%  of  the  expected  value;  well  within  acceptable  limits. 

//  Test  subroutine:  four_point 

std::cout  <<  "Test  subroutine:  f our_point\n\n" ; 
real  weight; 

std::cout  <<  "Test  #1:  All  4  points  in  Quadrant  I\n"; 
weight  = 

four_point (  0.001,0.004,  0.005,0.008,  0.007,0.001,  0.008,0.009,  0.0); 
std::cout  <<  "0.0002063887  "  <<  weight  <<  "\n\n"; 

std : : cout 

<<  "Test  #2:  2  points  in  Quadrant  I,  1  point  in  Quadrant  III\n"; 
weight  = 

four_point (-0.001, 0.004,  0.005,0.008,  0.007,0.001,  0.008,0.009,  0.0); 
std::cout  <<  "0.000253308876  "  <<  weight  <<  "\n\n"; 

std::cout  <<  "\n\n\n"; 

The  following  piece  of  code  was  written  for  nothing  more  than  to  be  able  to  view  the 
output  prior  to  the  main  test  program  completing.  Nothing  is  done  with  the  variable 
“something”. 

//  Pause  to  look  at  output  data 
int  something; 
std::cin  >>  something; 

return  0; 

}  //  main  subroutine  end 
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Subroutine:  power_avg 


The  subroutine  power_avg  integrates  the  area  under  the  curve  via  a  discrete 
trapezoidal  integration  scheme.  The  function  integrates  the  area  underneath  the  curve 
defined  in  the  power_funct  subroutine  from  time  t=tl  to  t=t2.  The  area  is  then  returned 
by  the  subroutine.  This  subroutine  was  not  used  in  the  final  program,  but  is  included  here 
for  follow-on  work  where  temporal  simulations  may  require  it  to  be  run. 

/* 

Subroutine:  power_avg 

Modified:  03/29/2005 

Created:  03/21/2005 

Creator:  Capt  Timothy  R.  Klein 

Description:  Integrates  the  area  under  the  curve  from  time  tl  to  time 
t2  . 

'k'k'k'k'k'k'k'k'ki<'k'k'ki<'k'k'ki<'k'k'k'k'k-k'ki<'k'k'k'k'k'k'ki<'k'k'ki<'k'k'ki<'k'k'k'k'k'k'ki<'k'k'ki<'k'k'ki<'k'k'k'k'k'k'k'k'k'k'ki<'k 

*/ 

real  power_avg (real  tl,real  t2) 

{ 

/*  Declare  Local  Variables  */ 
real  avg  =  0.0; 

real  num_steps  =  100.0; 

real  dt  =  1 . 0/ (num_steps*FREQ) ;  /*  Seconds  */ 

real  t; 

for (t=t 1 ; t<t2 ; t+=dt ) 

{ 

/*  Trapezoidal  Numerical  Integration  */ 

avg  +=  dt*(0.5*(  power_f unct (t ) +power_f unct (t+dt )  )); 

} 

return  avg; 

}  /*  power_avg  subroutine  end  */ 
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Subroutine:  power_funct 

The  subroutine  power_funct  returns  the  value  of  a  weighted  sine  function  given  a 
time  passed  to  it.  The  subroutine  uses  a  function  that  indicates  a  pull,  giving  amplitude  of 
- 1  for  the  first  part  of  the  wave.  The  second  part  of  the  wave  indicates  a  push,  giving 
amplitude  of  -7.  This  is  consistent  with  Font  [8],  where  there  is  7  times  more  push  than 
pull  during  a  cycle  containing  one  full  waveform.  This  subroutine  was  not  used  in  the 
final  program,  but  is  included  here  for  follow-on  work  where  temporal  simulations  may 


require  it  to  be  run. 


/* 

i(i<ici('ki<'k-k'k'k'k'k'k'k'k-k'ki<'k'k'k'k'k-k'ki<i('kic'kic'ki('kic'k'k'k'k-k'ki<'k'kic'kic'k'ki<ic'ki('kic'ki('k'k'k'k'kic'kici<ic'kic'kic 

Subroutine:  power_funct 


Modified: 
Created : 
Creator : 


03/29/2005 

03/21/2005 

Capt  Timothy  R.  Klein 


Description:  Power  delivery  function  represented  as  a  sine  wave  of 
varying  peaks . 

*/ 


real  power_funct (real  t) 

{ 

real  const  pi  =  3.14159;  /*  Pi,  constant 


V 


real  SINE_K1,  SINE_K2; 
real  w; 
real  t_temp; 
real  f_t  =  0.0; 


/*  Sine  wave  amplitudes  */ 

/*  Hz  frequency  and  Angular  frequency  */ 


/*  function  value  */ 

w  =  2.0*pi*FREQ;  /*  Angular  frequency  */ 

SINE_K1  =  -1;  /*  negative  first  half  of  sine  wave,  for  a  pull  */ 
SINE_K2  =  -7;  /*  positive  second  half  of  sine  wave,  for  a  push  */ 


t_temp  =  fmod (t, (1 . 0/FREQ) ) ; 
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/*  First  half  of  sine  wave  */ 

if (  t_temp>=0 . 0  &&  t_temp<= ( 0 . 5/FREQ)  ) 

{ 

f_t  =  (SINE_K1/ (SINE_K1+SINE_K2) ) *sin (w*t) ; 

} 

/*  Second  half  of  sine  wave  */ 

else  if (  t_temp> ( 0 . 5/FREQ)  &&  t_temp< ( 1 . 0/FREQ)  ) 

{ 

f_t  =  (SINE_K2/ (SINE_K1  +  SINE_K2)  ) *sin (w*t)  ; 

} 

return  f_t; 

}  /*  power_funct  subroutine  end  */ 
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Subroutine:  coord_xform 


The  subroutine  coord_xform  modifies  two  points  that  are  given  to  it,  pu  and  pv. 
This  is  done  by  first  finding  the  different  angles  between  the  x-axis  and:  the  point  to  be 
transformed,  the  transformation  vector,  and  finally  the  coordinates  of  the  tail  of  the 
transformation  vector.  Once  done,  all  of  the  points  and  vectors  are  rotated  by  the  angle 
that  was  found  for  the  transformation  vector.  The  transformed  point  has  the  transformed 


base  point  subtracted  from  it  as  it  is  now  the  new  (0,0)  point. 


/* 

Subroutine:  coord_xform 


Modified: 
Created : 
Creator : 


08/08/2005 

03/29/2005 

Capt  Timothy  R.  Klein 


Description:  Transforms  coordinates  and  vectors  to  give  2  sets  of 
points . 

ici<ici<ici<ici<ici<ici<ici<ici(ici<ici<ici<'k'kic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'kic'kic'k'k'k'k'k'k'k'k'kic'k'k'k'ki<'k'k'ki<'k'kic'k'k 

*/ 


void  coord_xf orm (real 

px. 

real 

py,  /*Coordinate  to  be  Transformed 

*/ 

real 

*pu. 

real 

*pv,  /‘Transformed  coordinate 

*/ 

real 

vlx , 

real 

vly,  /‘Vector  for  transformation 

*/ 

real 

plx. 

real 

ply  ) /‘Base  of  Vector 

*/ 

{ 

*pu  =  vlx* (px-plx)  +  vly* (py-ply) ; 
*pv  =  -vly* (px-plx)  +  vlx* (py-ply) ; 


return ; 

}  /*  coor_xform  subroutine  end  */ 
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Subroutine:  four_point 

The  subroutine  four_point  receives  4  points  in  no  particular  order  and  finds  which 
sets  of  2  points  are  adjacent  to  each  other  by  finding  the  longest  length  between  all  of  the 
points.  The  points  are  then  grouped  into  2  sets  of  3  points;  each  set  shares  2  points  that 


were  adjacent  to  each  other.  These  points  are  sent  to  the  subroutine  three-point  for  line- 


straddle  checking  and  finally  integration. 


/* 

Subroutine:  four_point 


Modified: 
Created : 
Creator : 


08/08/2005 

03/21/2005 

Capt  Timothy  R.  Klein 


Description:  Takes  in  4  points  in  no  particular  order  and  finds 

adjacent  points.  The  points  are  then  grouped  into  3's, 
sharing  2  opposite  sets  of  points,  and  sent  to  the 
three_point  subroutine  for  integration. 

*/ 

real  four_point (  real  ul, real  vl, 
real  u2,real  v 2, 
real  u3,real  v3, 
real  u4,real  v4, 
real  weight) 


/*  Declare  Local  Point  Variables  */ 
real  temp_weight  =  0.0; 


/*  Calculate  the  vector  magnitudes  for  all  6  vectors  from  4  points  */ 


real 

vl2_mag  =  sqrt ( 

pow ( (u2-ul ) 

,2) 

+ 

pow ( ( v2-vl 

,2) 

)  ; 

real 

vl3_mag  =  sqrt ( 

pow ( (u3-ul ) 

,2) 

+ 

pow ( ( v3-vl 

,2) 

)  ; 

real 

vl4_mag  =  sqrt ( 

pow ( (u4-ul ) 

,2) 

+ 

pow ( ( v4-vl 

,2) 

)  ; 

real 

v23_mag  =  sqrt ( 

pow ( (u3-u2 ) 

,2) 

+ 

pow ( ( v3-v2 

,2) 

)  ; 

real 

v24_mag  =  sqrt ( 

pow ( (u4-u2 ) 

,2) 

+ 

pow ( ( v4-v2 

,2) 

)  ; 

real 

v34_mag  =  sqrt ( 

pow ( (u4-u3 ) 

,2) 

+ 

pow ( ( v4-v3 ' 

,2) 

) ; 
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/*  Determine  the  maximum  magnitude  */ 
real  maximum  =  vl2_mag; 

if (  maximum  <  vl3_mag  )  maximum=vl3_mag; 
if (  maximum  <  vl4_mag  )  maximum=vl 4_mag; 
if (  maximum  <  v23_mag  )  maximum=v23_mag; 
if (  maximum  <  v24_mag  )  maximum=v2 4_mag; 
if (  maximum  <  v34_mag  )  maximum=v34_mag; 


/*  Determine  opposite  points  */ 

if (  maximum==vl2_mag  | |  maximum==v34_mag  ) 

{ 

temp_weight  +=  three_point (  ul,vl,  u3,v3,  u4,v4,  temp_weight  ) 
temp_weight  +=  three_point (  u2,v2,  u3,v3,  u4,v4,  temp_weight  ) 

} 

else  if (  maximum==vl3_mag  | |  maximum==v2 4_mag  ) 


temp_weight  +=  three_point (  ul,vl,  u2,v2,  u4,v4, 
temp_weight  +=  three_point (  u3,v3,  u2,v2,  u4,v4, 

} 

else  if (  maximum==vl 4_mag  | |  maximum==v23_mag  ) 

{ 

temp_weight  +=  three_point (  ul,vl,  u2,v2,  u3,v3, 
temp_weight  +=  three_point (  u4,v4,  u2,v2,  u3,v3, 

} 


temp_weight  ) 
temp_weight  ) 


temp_weight  ) 
temp_weight  ) 


return  temp_weight; 

}  /*  four_point  subroutine  end  */ 
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Subroutine:  three_point 

The  subroutine  three_point  receives  3  points  in  no  particular  order  and  checks  if 
the  area  enclosed  by  the  points  straddles  any  boundary  lines  set  by  the  user.  If  the  points 
are  all  contained  within  a  boundary  region,  then  volume  integration  is  completed  and  the 
weight  is  returned. 

Otherwise,  the  line  and  boundary  intercepts  are  found.  A  new  set  of  3  points  is 
generated  along  with  a  set  of  4  points.  Each  set  should  lie  on  opposite  sides  of  the 


boundary.  The  subroutine  is  recursive  until  all  boundaries  have  not  been  crossed  by  a  set 
of  points. 


/* 

Subroutine:  three_point 


Modified: 
Created : 
Creator : 


08/08/2005 

03/21/2005 

Capt  Timothy  R.  Klein 


Description:  Takes  in  3  points  in  no  particular  order  and  looks  to  see 
if  the  area  enclosed  by  the  points  contains  any  boundaries 
defined  in  x_cross  and  y_cross.  If  there  is  a  boundary, 
then  the  area  will  be  parsed  into  a  3  point  and  4  point  set 
to  undergo  evaluation  again.  If  there  is  no  boundary,  then 
the  3  points  containing  the  area  will  be  passed  to  a 
numerical  integrator  to  return  the  weight  or  volume  in  that 
area . 


*/ 

real  three_point (real  ul, real  vl, 
real  u2,real  v 2, 
real  u3,real  v3, 
real  weight) 


/*  Declare  Local  Point  Variables  */ 
int  check  =  0; 

real  alu,  alv,  a2u,  a2v;  /  *  Points  on  line  *  / 
real  temp_weight  =  0.0; 

int  n; 
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for  (n=l; n<6; ++n) 


{ 

/*  Check  to  see  if  points  are  on  both  sides  of  line  */ 
if (  (  line_side (n, ul , vl ) <0 . 0  && 

line_side (n, u2 , v2 ) <0 . 0  && 
line_side (n, u3 , v3 ) >0 . 0  )  | | 

(  line_side (n,  ul ,  vl ) <0 . 0  && 

line_side (n, u2 , v2 ) >0 . 0  && 
line_side (n, u3 , v3 ) <0 . 0  )  | | 

(  line_side (n,  ul ,  vl ) >0 . 0  && 

line_side (n, u2 , v2 ) <0 . 0  && 
line_side (n, u3 , v3 ) <0 . 0  )  |  | 


(  line_side (n, ul , vl ) >0 . 0  && 
line_side (n, u2 , v2 ) >0 . 0  && 
line_side (n, u3 , v3 ) <0 . 0  )  | | 

(  line_side (n, ul ,  vl ) >0 . 0  && 
line_side (n, u2 , v2 ) <0 . 0  && 
line_side (n, u3 , v3 ) >0 . 0  )  |  | 

(  line_side (n, ul , vl ) <0 . 0  && 
line_side (n, u2 , v2 ) >0 . 0  && 
line_side (n, u3 , v3 ) >0 . 0  )  ) 


check  =  1;  /*  area  being  parsed,  no  need  to  run  integration  */ 

/*  pi  and  p2  on  same  side  of  line  */ 

if (  line_side (n, ul , vl ) *line_side (n, u2 , v2 )  >  0.0  ) 

{ 

/*  Using  p3  as  the  vertex,  */ 

/*  find  (u, v) -coordinates  where  lines  intercept  */ 
line_intercept (n,  u3,v3,  ul,vl,  &alu,&alv); 
line_intercept (n,  u3,v3,  u2,v2,  &a2u,&a2v); 


temp_weight  +=  four_point  (  ul,  vl, 

alu, alv, 

temp_weight  +=  three_point (  u3,  v3, 

alu, alv. 


} 


u2,  v2, 
a2u, a2v, 

a2u, a2v. 


temp_weight  ) ; 
temp_weight  ) ; 


/*  pi  and  p3  on  same  side  of  line  */ 

else  if (  line_side (n,ul, vl) *line_side (n,u3,  v3)  >  0.0  ) 

{ 

/*  Using  p2  as  the  vertex,  */ 

/*  find  (u, v) -coordinates  where  lines  intercept  */ 
line_intercept (n,  u2,v2,  ul,vl,  &alu,&alv); 
line_intercept (n,  u2,v2,  u3,v3,  &a2u,&a2v); 


temp_weight  +=  four_point  (  ul,  vl, 

alu, alv, 

temp_weight  +=  three_point (  u2,  v2, 

alu, alv. 


} 


u3,  v3, 

a2u,a2v,  temp_weight  ); 
a2u, a2v,  temp_weight  ) ; 
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/*  p2  and  p3  on  same  side  of  line  */ 

else  if (  line_side (n, u2, v2) *line_side (n,  u3,  v3)  >  0.0  ) 

{ 

/*  Using  pi  as  the  vertex,  */ 

/*  find  (u, v) -coordinates  where  lines  intercept  */ 

1 ine_intercept (n,  ul,vl,  u2,v2,  &alu,&alv); 
line_intercept (n,  ul,vl,  u3,v3,  &a2u,&a2v); 


temp_weight 

temp_weight 


} 


+=  four_point  (  u2,  v2, 
alu, alv, 
+=  three_point (  ul,  vl, 
alu, alv. 


u3,  v3, 

a2u,a2v,  temp_weight  ); 
a2u,a2v,  temp_weight  ); 


}  /*  line_side  if  check  */ 

}  /*  n  for  loop  */ 

/*  Area  was  not  split  along  any  boundaries,  */ 

/*  integrate  area  for  volume  weight  */ 

if (check  ==  0) 

{ 

/*  Passed  all  checks  on  boundaries,  */ 

/*  can  now  integrate  area  WRT  function  */ 

temp_weight  +=  volume_integration (ul , vl ,  u2,v2,  u3,v3,  temp_weight) 

} 


return  temp_weight; 

}  /*  three_point  subroutine  end  */ 
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Subroutine:  line_side 


The  subroutine  line_side  receives  a  point  as  well  as  the  line  number  to  compare 
against.  The  line  locations  are  defined  in  the  subroutine  line_offset.  The  subroutine 
line_side  outputs  the  distance  from  the  line  in  either  x  or  y  coordinates,  depending  upon 


the  line  direction.  If  the  number  to  be  returned  is  negative,  then  the  point  is  to  either  the 


left  or  below  the  line.  The  opposite  is  true  for  a  number  that  is  returned  positive. 


/* 

Subroutine:  line_side 


Modified: 
Created : 
Creator : 


08/08/2005 

03/28/2005 

Capt  Timothy  R.  Klein 


Description:  Gives  the  side  of  the  line  the  point  (u,v)  is  on. 

(-)  is  the  left  side,  (+)  is  the  right  side 

'k'k'k'k'ki<'k'k'k'k'k'k'ki<'k'k'k'k'k'k'ki<'k'k'ki<'k'k'ki<'k'k'k'k'k'k'k'k'k-k'ki<'k'k'ki<'k'k'ki<'k'k'k'k'k'k'ki<'k'k'k'k'k'k'k'k'k-k'kj<'k 

*/ 

real  line_side (int  n,  real  u,  real  v) 


real  offset; 
real  side  =  0.0; 


offset  =  line_offset (n) ; 

switch (n) 

{ 

case  1  : 
case  2  : 
case  3  : 
case  4  : 

side  =  u  -  offset; 
break ; 
case  5  : 
case  6: 

side  =  v  -  offset; 
break ; 

} 

return  side; 

}  /*  line_side  subroutine  end  */ 
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Subroutine:  line_intercept 

The  subroutine  line_intercept  receives  a  set  of  points  as  well  as  the  line  number  to 
compare  against.  Because  this  function  is  set  as  void,  there  is  no  return  value.  Instead, 
the  changes  made  by  the  function  are  done  directly  to  the  pointers  au  and  av.  The  line 
locations  are  defined  in  the  subroutine  line_offset.  The  points  au  and  av  are  the 
coordinates  of  where  the  intersection  of  the  line  formed  by  (ul,vl)  and  (u2,v2)  and  the 


line  defined  by  n  when  passed  to  the  subroutine  line_offset. 


/* 

Subroutine:  line_intercept 


Modified: 
Created : 
Creator : 


08/08/2005 

03/28/2005 

Capt  Timothy  R.  Klein 


Description:  Gives  the  cross-product  of  a  vector  (vx,vy)  with  tail  at 
point  (px,py)  and  the  resulting  vector  (x-px,y-py)  when 
point  (x,y)  is  the  head  and  point  (px,py)  is  the  tail.  For 
simplicity,  (px, py) = (0, 0) ,  thus  making  the  resulting  vector 
(x,  y) 


void  line_intercept (int  n 

r 

/* 

line  number 

*/ 

real 

ul, 

real 

vl. 

/* 

point  pi  =  (ul,vl) 

*/ 

real 

u2 , 

real 

v2 , 

/* 

point  p2  =  (u2,v2) 

*/ 

real 

*au, 

real 

*av) 

/* 

intercept  point  (au, av) 

*/ 

{ 

real  offset; 

real  ul_new, u2_new,  vl_new, v2_new; 


offset  =  line_offset (n)  ; 


switch (n) 

{ 

case  1  : 
case  2  : 
case  3  : 
case  4  : 

ul_new  =  ul  -  offset;  /*  location  WRT  line  */ 
u2_new  =  u2  -  offset;  /*  location  WRT  line  */ 
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if (  (ul_new  <=  0.0  &&  u2_new  >=0.0)  | | 

(ul_new  >=  0.0  &&  u2_new  <=  0.0)  ) 

{ 

*au  =  0.0  +  offset; 

*av  =  vl  +  (ul_new/ (ul_new-u2_new) ) * (v2-vl )  ; 

} 

else 

{ 

*au  =  0.0; 

*  av  =  0.0; 

} 

break; 


case  5  : 
case  6: 

vl_new  =  vl  -  offset;  /’ 
v2_new  =  v2  -  offset;  /’ 
if (  (vl_new  <=  0.0  &&  v2_new  >=0.0)  || 
(vl_new  >=  0.0  &&  v2_new  <=  0.0)  ) 


location  WRT  line  */ 
location  WRT  line  */ 


*av  =  0.0  +  offset; 

*au  =  ul  +  (vl_new/ (vl_new-v2_new) ) * (u2-ul ) ; 

} 

else 

{ 

*au  =  0.0; 

*av  =  0.0; 

} 

break; 


return ; 

/*  line_intercept  subroutine  end  */ 
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Subroutine:  line_offset 


Line  Number  Reference  Figure 


The  subroutine  line_offset  receives  the  line  number  and  returns  the  value  of  the 


location  of  that  line.  This  is  an  easy  place  to  modify  a  globally  used  variable  when 
making  modifications.  Parameterizations  of  experiments  show  that  the  DBD  main  wall 
jet  effects  are  seen  within  1  cm  of  the  DBD  device,  with  a  glow  surrounding  the 


electrodes  as  well.  Thus,  the  following  bounding  for  the  wall  jet  was  assumed. 


/* 

Subroutine:  line_offset 


Modified: 
Created : 
Creator : 


08/09/2005 

03/28/2005 

Capt  Timothy  R.  Klein 


Description:  Quick  set  area  to  line  offsets. 
*/ 

real  line_of f set ( int  n) 


real  offset  =  0.0; 

switch (n) 

{ 
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case  1:  /*  linelx  */ 

offset  =  -0.002;  /*  Original 
switch (WScheme) 

{ 


case  1  : 
offset 
break; 
case  2  : 
offset 
break; 
default : 
offset 
break; 


/*  Bouef f SPitchf ord 
=  -0.0001;  /*  -100 

/*  Roy&Gaitonde 
=  -0 .0135;  /*  -1 . 35 

/*  Original 
=  -0.0020; 


break; 


case  2:  /*  line2x  */ 

offset  =  0.000;  /*  Original 

switch (WScheme) 

{ 


case  1  : 
offset 
break; 
case  2  : 
offset 
break; 
default : 
offset 
break; 


/*  Bouef f SPitchf ord 

=  0.0000; 

/*  Roy&Gaitonde 
=  -0.0045; 

/*  Original 

=  0.0000; 


break ; 


case  3:  /*  line3x  */ 

offset  =  +0.010;  /*  Original 
switch (WScheme) 

{ 


case  1  : 
offset 
break; 
case  2  : 
offset 
break; 
default : 
offset 
break; 

} 

break ; 


/*  Bouef f SPitchf ord 
=  +0.0006;  /*  600 

/*  Roy&Gaitonde 
=  +0.0040;  /*  0.4 

/*  Original 

=  +0 .010;  /*  1.0 


*/ 

*/ 

urn  */ 
*/ 

cm  */ 
*/ 


*/ 

*/ 

*/ 

*/ 


*/ 

*/ 

urn  *  / 
*/ 

cm  *  / 
*/ 

cm  */ 
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case  4:  /*  line4x  */ 

offset  =  +0.015;  /*  Original  */ 
switch (WScheme) 

{ 


case  1  : 
offset 
break; 
case  2  : 
offset 
break; 
default : 
offset 
break; 


/*  Bouef f &Pitchf ord  */ 

=  +0.0006;  /*  600  um  */ 

/*  Roy&Gaitonde  */ 

=  +0 .0165;  /*  1.65  cm  */ 

/*  Original  */ 

=+0.010;  /*  1.00  cm  */ 


break; 


case  5:  /*  linely  */ 

offset  =  0.000;  /*  Original  */ 

switch (WScheme) 

{ 


case  1  : 

/* 

Bouef f&Pitchf ord 

*/ 

offset 

break; 

= 

0 . 000; 

case  2  : 

/* 

Roy&Gaitonde 

*/ 

offset 

break; 

= 

0 . 000; 

default : 

/* 

Original 

*/ 

offset 

break; 

= 

0 . 000; 

} 

break ; 


case  6:  /*  line2y  */ 

offset  =  +0.003;  /*  Original 
switch (WScheme) 

{ 


case  1  : 
offset 
break; 
case  2  : 
offset 
break; 
default : 
offset 
break; 

} 

break ; 


/*  Bouef f&Pitchf ord 
=  +0.00005;  /*  50 

/*  Roy&Gaitonde 
=  +0.005;  /*  0.5 

/*  Original 
=  +0.003;  /*  0.3 


*/ 

*/ 

um 

*/ 

cm 

*/ 

cm 


*/ 

*/ 


*/ 


return  offset; 

}  /*  line_offset  subroutine  end  */ 
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Subroutine:  curve_y 

The  subroutine  curve_y  receives  the  x- coordinate  for  a  point  on  a  curve  of  the 
airfoil.  The  curve  may  be  flat,  as  in  the  case  of  the  flat-plate.  However  it  may  have 
curvature  as  described  by  this  subroutine.  The  variable  cord_length  is  used  as  a  scaling 
factor  for  the  location  of  the  x-coordinate. 


/* 

Subroutine:  curve_y 


Modified: 
Created : 
Creator : 


08/08/2005 

03/21/2005 

Capt  Timothy  R.  Klein 


Description:  Given  a  symmetric  airfoil  and  a  x-coordinate,  function 

will  return  y.  cord_length  is  used  as  a  scaling  factor. 


*/ 

real  curve_y (real  x_in,  real  cord_length) 


real  y; 


real  const  D_1  =  0.0; 

real  const  D_2  =  0.0; 

real  const  D_3  =  0.0; 

real  const  D_4  =  9.0; 

if (F_OR_A  ==  0)  /*  Flat-Plate  */ 

{ 

real  offset  =  0 . 5* (cord_length-CORD)  -  Y_POSITION*CORD; 
y  =  0 . 0  +  offset ; 

} 

else  if(F_OR_A  ==  1)  /*  Airfoil*/ 

{ 

real  f  =  D_1 ;  /*  Maximum  Camber  */ 

real  xf  =  D_2 ;  /*  Position  of  Maximum  Camber  */ 

real  t  =  0.1*D_3  +  0.01*D_4;  /*  %  thickness/cord  */ 

/*  return  to  normalized  x-coordinates*/ 
real  x  =  x_in  /  cord_length; 

real  c  =  1.0;  /*  Normalized  Cord  Length  */ 

real  xl  =  xf/c; 
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real  yc  =  c* 

(f/c) * 

pow (1 . 0/ (1 . 0-xl) , 2) * 

(  (1.0-2.0*xl)  +  2.0*xl*(x/c)  -  (x*x) / (c 

real  yt  =  c*5*t*(  +  (  0 . 29690*pow (x, 0 . 5)  ) 

-  (  0 . 12600*pow (x, 1  )  ) 

-  (  0 . 35160*pow (x, 2  )  ) 

+  (  0 . 28430*pow  (x, 3  )  ) 

-  (  0 . 10150*pow (x, 4  )  )  ); 

/*  y-coordinate  with  cord_length  scaling  factor  */ 
y  =  cord_length* (yc  +  yt); 


return  y; 

/*  curve_y  subroutine  end  */ 
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Subroutine:  curve_dy 


The  subroutine  curve_dy  receives  the  x- coordinate  for  a  point  on  a  curve  of  the 
airfoil.  The  curve  may  be  flat,  as  in  the  case  of  the  flat-plate.  However  it  may  have 
curvature  as  described  by  this  subroutine.  Either  case  will  change  the  output  given  the  x- 


coordinate,  as  this  subroutine  returns  the  slope  at  this  point.  The  variable  cord_length  is 


used  as  a  scaling  factor  for  the  location  of  the  x- coordinate. 


/* 

Subroutine:  curve_dy 


Modified: 
Created : 
Creator : 


04/05/2005 

03/28/2005 

Capt  Timothy  R.  Klein 


Description:  Given  a  symmetric  airfoil  and  a  x-coordinate,  function 
will  return  slope= (dy/dx) .  cord_length  is  used  as  a 
scaling  factor. 


*/ 

real  curve_dy (real  x_in,  real  cord_length) 


real  dy_dx; 


real  const  D_1  =  0.0; 
real  const  D_2  =  0.0; 
real  const  D_3  =  0.0; 
real  const  D_4  =  9.0; 


if (F_OR_A  ==  0)  /*  Flat-Plate  */ 


dy_dx  =  0.0; 

} 

else  if (F_OR_A  ==  1)  /*  Airfoil  */ 

{ 

real  f  =  D_1 ;  /*  Maximum  Camber  */ 

real  xf  =  D_2 ;  /*  Position  of  Maximum  Camber  */ 

real  t  =  0.1*D_3  +  0.01*D_4;  /*  %  thickness/cord  */ 


/*  return  to  normalized  x-coordinates  */ 
real  x  =  x_in  /  cord_length; 

real  c  =  1.0;  /*  Normalized  Cord  Length  */ 


real  xl  =  xf/c; 
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real  dyc_dx  =  c* 

(f/c) * 

pow  ( 1 . 0 / ( 1 . 0-xl ) , 2 ) * 

(  (1.0-2. 0  *xl )  +  2 . 0  *xl *  (1.0/c)  -  (2 . 0*x) / (c*c)  ); 

real  dyt_dx  =  c*5*t*(  +  (  0 . 5*0 . 29690*pow (x, -0 . 5)  ) 

-  (  0.12600  ) 

-  (  2 . 0*0 . 35160*pow  (x,  1  )  ) 

+  (  3 . 0*0 . 28430*pow (x,  2  )  ) 

-  (  4 . 0*0 . 10150*pow (x,  3  )  )  ); 

/*  slope  of  line  given  x-coordinate  */ 

/*  normalized  with  cord_length  scaling  factor  */ 
dy_dx  =  (dyc_dx  +  dyt_dx) ; 

} 

return  dy_dx; 

}  /*  curve_dy  subroutine  end  */ 
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Subroutine:  volume_integration 

The  subroutine  volume_integration  receives  3  points  that  do  not  cross  over  any 


boundaries  and  returns  the  weight  via  discrete  integration. 


/* 

■k'k'kit'kit'kit'kit'kit'kit'kir'kit'kit'kit'kit'k'k'kit'k'k'kit'k'k'kit'k'k'kit'k'k'kit'k'k'kit'k'k'kit'k'k'kit'kit'kit'k'k'kit'k'k'kit'kit-k 

Subroutine:  volume_integration 


Modified: 
Created : 
Creator : 


04/08/2005 

03/22/2005 

Capt  Timothy  R.  Klein 


Description:  Integrates  the  area  under  the  curve  for  the  area  contained 
within  the  3  given  points. 

*/ 

real  volume_integration (real  ul, real  vl, 

real  u2,real  v2, 
real  u3,real  v3, 
real  weight) 


/*  Declare  Local  Variables  */ 
real  s,ds; 
real  t,dt; 

real  plu,plv,  p2u,p2v,  p3u,p3v,  p4u,p4v; 
real  temp_weight  =  0.0; 


int  i, j; 

/*  Number  of  divisions  in  cell  area  =  0 . 5*num_div* (num_div+l )  */ 

int  num_div  =  50; 


/*  Vector  1  */ 
real  vlu  =  u2-ul; 
real  vlv  =  v2-vl; 


/*  Vector  2  */ 
real  v2u  =  u3— ul; 
real  v2v  =  v3-vl; 


real  N  =  abs (vlu*v2v  -  v2u*vlv) ; 
s  =  0.0; 

for ( i=l ; i<=num_div; ++i ) 

{ 


ds  =  ( 1 . 0  /  (real )  (num_div) ) * ( 1 . 0-0.0); 
t  =  0.0; 
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f or  (  j=l; j<= (num_div-i  +  l )  ;  ++ j  ) 

{ 


dt  =  ( 1 . 0/ ( real )  (num_div-i  +  l ) ) * ( 1 . 0-s )  ; 


if (  (t+dt ) < ( 1 . 0-s )  )  /*  Square  area's  in  the  middle  */ 

{ 

/*  Coordinate  Transformations  */ 


plu 

=  (vlu*  (s  ) 

+ 

v2u*  (t 

)  ) 

+ 

ul; 

plv 

=  (vlv* (s  ) 

+ 

v2v* (t 

)  ) 

+ 

vl; 

p2u 

=  (vlu* ( s+ds ) 

+ 

v2u*  (t 

)  ) 

+ 

ul; 

p2v 

=  (vlv* (s+ds ) 

+ 

v2v* (t 

)  ) 

+ 

vl; 

p3u 

=  (vlu*  (s  ) 

+ 

v2u* (t+dt ) ) 

+ 

ul; 

p3v 

=  (vlv* (s  ) 

+ 

v2v* (t+dt ) ) 

+ 

vl; 

p4u 

=  (vlu* (s+ds ) 

+ 

v2u* (t+dt ) ) 

+ 

ul; 

p4v 

=  (vlv* (s+ds ) 

+ 

v2v* (t+dt ) ) 

+ 

vl; 

temp 

_weight  +=  N*(l. 

. 0*ds*dt) 

*0  . 

25 

* 

(  weight_f unct (plu, plv)  + 
weight_funct (p2u, p2v)  + 
weight_funct (p3u, p3v)  + 
weight_funct (p4u, p4v)  ); 

} 

/*  Triangle  area's  along  the  s=l-t  edge  */ 
else  if (  (t+dt )>=( 1 . 0-s )  ) 


/*  Coordinate  Transformations  */ 


plu 

plv 

p2u 

p2v 

p3u 

p3v 


)  +  v2u* (t 

)  +  v2v* (t 


(vlu*  (s 
(vlv* ( s 
(vlu* (s+ds)  +  v2u*  (t 
(vlv* (s+ds)  +  v2v* (t 


) )  +  ul ; 
) )  +  vl ; 
) )  +  ul ; 
) )  +  vl ; 


(vlu*  (s  )  +  v2u*(t+dt))  +  ul; 
(vlv* (s  )  +  v2v*(t+dt))  +  vl ; 


temp_weight  +=  N* (0 . 5*ds*dt) * (1 . 0/3 . 0) * 

(  weight_funct (plu,  plv)  + 
weight_funct (p2u, p2v)  + 
weight_funct (p3u, p3v)  ); 

} 

t  +=  dt ; 

}  /*  j  loop  */ 

s  +=  ds; 

}  /*  i  loop  */ 

return  temp_weight; 

}  /*  volume_integration  subroutine  end  */ 


VH-42 


Subroutine:  weight_funct 

The  subroutine  weight_funct  receives  a  point  and  returns  its  value.  This 
subroutine  was  necessary  for  fast  modification  of  the  weighting  function  equation  and 
how  it  acted  on  all  sides  of  the  boundaries. 


/* 

'ki<'k'k'k'k'k'k'k'k'k'k'ki<'k'k'k'k'k'k'ki<'k'k'k'k'k'k'ki<'k'k'k'k'k'k'k'k'k'k'k'k'k'k'ki<'k'k'ki<'k'k'ki<'k'k'ki<'k'k'ki<'k'k'k'k'k'k'k'k'k 

Subroutine:  weight_funct 


Modified: 
Created : 
Creator : 


08/09/2005 

03/22/2005 

Capt  Timothy  R.  Klein 


Description:  2D  function  of  cell  weight  dependent  on  position. 
*/ 

real  weight_f unct ( real  u,  real  v) 


/*  Declare  Local  Variables  */ 

real  f_xy;  /*  function  value 

real  e,p,mid; 


real  K_l, 
real  K_2, 
real  K_3, 


C_1U,  C_1V,  C_1 ; 
C_2U,  C_2V,  C_2 ; 
C_3U,  C_3V,  C_3 ; 


*/ 


/* 

Switch  for  Bouef f &Pitchf ord  or  Roy&Gaitonde  Weighting-Scheme 
Bouef f &Pitchford:  WScheme  =  1 
Roy&Gaitonde:  WScheme  =  2 

*/ 


switch (WScheme) 

{ 
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case  1:  /*  Bouef f SPitchf ord  */ 

/*  First  Function  Weighting  Constants  */ 

K_1  =  +  10000.0;  /*  Weight  of  Amplitude  */ 

C_1U  =  +  50000.0;  /*  Decent  Rate  constants  */ 

C_1V  =  +  80000.0;  /*  Decent  Rate  constants  */ 

/*  Second  Function  Weighting  Constants  */ 

K_2  =  +  10000.0;  /*  Weight  of  Amplitude  */ 

C_2U  =  +  100.0;  /*  Decent  Rate  constants  */ 

C_2V  =  +  80000.0;  /*  Decent  Rate  constants  */ 


/*  -0.0001  <=  u  <  0.000  */ 

if (  u  >=  line_of f set ( 1 )  &&  u  <  line_of f set (2)  ) 


/*  -0.001  <=  v  <=  +0.00005  */ 

if (  v  <=  line_of f set (6)  &&  v  >=  -0.001) 

{ 

if (u<0 . 0 )  u=-u; 
if (v<0 . 0)  v=-v; 
e  =  0.0  -  u*C_lU  -  v*C_lV; 
f_xy  =  (exp (e ) ) *K_1 ; 

} 

else 


f_xy  =  0.0; 

} 

} 


/*  0.000  <=  u  <=  +0.00060  */ 

else  if (  u  >=  line_of f set (2)  &&  u  <=  line_offset (3)  ) 


/*  -0.001  <=  v  <=  +0.00005  */ 
if (  v  <=  line_of f set (6)  &&  v  >=  -0.001) 

{ 

if (u<0 . 0 )  u=-u; 
if (v<0 . 0 )  v=-v; 

e  =  0.0  -  u*C_2U  -  v*C_2V ; 
f_xy  =  (exp (e ) ) *K_2 ; 

} 

else 

{ 

f_xy  =  0.0; 

} 


else  /*  everywhere  else  f_xy=0 . 0  */ 

{ 

f_xy  =  0.0; 

} 

break ; 
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case  2:  /*  Roy&Gaitonde  */ 

/*  First  Function  Weighting  Constants  */ 

K_1  =  163.0  ;  /*  Weight  of  Amplitude  */ 

C_1  =  5.0  ;  /*  Decent  Rate  constants  */ 

/*  Second  Function  Weighting  Constants  */ 

K_2  =  3200.0  ;  /*  Weight  of  Amplitude  */ 

C_2  =  2.75;  /*  Decent  Rate  constants  */ 

/*  Third  Function  Weighting  Constants  */ 

K_3  =  100.0  ;  /*  Weight  of  Amplitude  */ 

C_3U  =  1.20*2;  /*  Decent  Rate  constants  */ 

C_3V  =  2.50*2;  /*  Decent  Rate  constants  */ 

/*  -0.0135  <=  u  <  -0.0045  */ 

if (  u  >=  line_of f set ( 1 )  &&  u  <  line_offset (2)  ) 

{ 

/*  -0.0010  <=  v  <=  +0.0050  */ 

if (  v  <=  line_offset (6)  &&  v  >=  -0.001) 

{ 

mid  =  (  line_of f set ( 1 )  +  line_of f set (2)  )  /  2.0; 

u  =  u  -  mid; 

p  =  0.0  -  C_l*sqrt (u*u+v*v) ; 

if (u<0 . 0 )  u=-u; 
if (v<0 . 0)  v=-v; 
p  =  0.0  -  C_l*u  -  C_l*v; 
f_xy  =  K_1 * (exp (p) ) ; 

} 

else 

{ 

f_xy  =  0.0; 

} 

} 

/*  -0.0045  <=  u  <=  +0.0040  */ 

else  if (  u  >=  line_of f set (2)  &&  u  <=  line_of f set (3)  ) 

{ 

/*  -0.0010  <=  v  <=  +0.0050  */ 

if (  v  <=  line_offset (6)  &&  v  >=  -0.001) 

{ 

mid  =  (  line_offset (2)  +  line_offset (3)  )  /  2.0; 

u  -=  mid; 

p  =  0.0  -  C_2*sqrt (u*u+v*v)  ; 

if (u<0 . 0 )  u=-u; 
if (v<0 . 0)  v=-v; 

p  =  0.0  -  u*100 . 0*C_2  -  v* 1 0 0 . 0 *C_2 ; 
f_xy  =  K_2*pow ( 10 . 0 , p) ; 

} 

else 

{ 

f_xy  =  0.0; 

} 
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/*  +0.0040  <  u  <=  +0.0165  */ 

else  if (  u  >  line_offset (3)  &&  u  <=  line_of f set (4)  ) 

{ 

/*  -0.0010  <=  V  <=  +0.0050  */ 

if (  v  <=  line_of f set (6)  &&  v  >=  -0.001) 


mid  =  (  line_offset (3)  +  line_of f set (4 )  ) 


u  -=  mid; 
if  (u<0 . 0)  u=-u; 
if (v<0 . 0 )  v=— v; 

p  =  0.0  -  u*100 . 0*C_3U  -  v*100  . 0*C_3V; 
f_xy  =  0.0  -  K_3 *pow ( 10 . 0 , p) ; 

} 

else 


/  2.0; 


f_xy  =  0.0; 

} 

} 

else  /*  everywhere  else  f_xy=0 . 0  */ 

{ 

f_xy  =  0.0; 

} 

break; 
default : 

printf ( "Weighting-Scheme  is  neither  Bouef f SPitchf ord  or 
Roy&Gaitonde " ) ; 


/*  First  Function  Weighting  Constants  */ 


K_ 

.1  = 

o 

\ — 1 

/* 

Weight 

of  Amplitude 

*/ 

C_ 

1U  = 

800.0; 

/* 

Decent 

Rate  constants 

*/ 

C_ 

.lv  = 

800 .0; 

/* 

Decent 

Rate  constants 

*/ 

/* 

Second  Function 

Weighting  Constants  */ 

K_ 

2  = 

7.0; 

/* 

Weight 

of  Amplitude 

*/ 

C_ 

2U  = 

100 . 0; 

/* 

Decent 

Rate  constants 

*/ 

C_ 

2  V  = 

800.0; 

/* 

Decent 

Rate  constants 

*/ 

/*  -0.002  <=  u  <  0.000  */ 

if (  u  >=  line_offset (2)  &&  u  <  line_of f set ( 1 )  ) 

{ 

/*  -0.001  <=  V  <=  0.010  */ 

if  (  v  <=  line_of f set (6)  &&  v  >=  -0.001) 

{ 

if (u<0 . 0 )  u=-u; 
if (v<0 . 0 )  v=-v; 
e  =  0.0  -  u*C_lU  -  v*C_lV; 
f_xy  =  (exp (e ) ) *K_1 ; 

} 

else 

{ 

f_xy  =  0.0; 


} 
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/*  0.000  <=  u  <=  0.010  */ 

else  if (  u  >=  line_of f set (1)  &&  u  <=  line_offset (3)  ) 
{ 

/*  -0.001  <=  v  <=  0.010  */ 

if (  v  <=  line_of f set (6)  &&  v  >=  -0.001) 

{ 

if (u<0 . 0 )  u=-u; 
if (v<0 . 0 )  v=-v; 

e  =  0.0  -  u*C_2U  -  v* C_2V; 
f_xy  =  (exp (e ) ) *K_2 ; 

} 

else 


f_xy  =  0.0; 

} 

} 

else  /*  everywhere  else  f_xy=0 . 0  */ 

{ 

f_xy  =  0.0; 

} 

break; 


return  f_xy; 

/*  weight_funct  subroutine  end  */ 
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VIII.  Appendix  C 


The  large  number  of  data  sets  and  resulting  simulations  to  each  set  required  an 
incredible  amount  of  processing.  To  alleviate  the  monotony  and  presence  of  the  operator 
in  the  computer  lab,  a  PBS  (Portable  Batch  System)  script  and  a  Journal  script  were 
created  for  each  simulation  to  automate  its  setup  and  execution.  This  allowed  for  a  quick 
initiation  of  the  simulation  from  a  remote  site. 

PBS  Script 

The  operator  first  remotely  logged- on  to  the  AFIT  Tahoe  cluster  computer  and 
moved  to  the  appropriate  simulation  set  directory.  A  check  by  typing  “qstaf  ’  was 
initiated  to  see  the  usage  of  the  machine  and  to  check  on  existing  jobs  already  submitted 
to  the  queue.  Submitting  a  job  to  the  queue  was  done  by  typing  “qsub”  and  the  PBS 
script  file  name.  This  was  then  executed  by  the  Tahoe  cluster  and  assigned  by  the  cluster 
to  a  node  for  processing.  The  following  PBS  script  is  the  template  example  that  was  used 
for  each  simulation. 

# ! /bin/bash 

#PBS  -o  out 

#PBS  -e  error 

#PBS  -1  nodes=l : ppn=2 

#PBS  -j  oe 

Cd  $PBS_0_W0RKDIR 

rm  -R  -f  libudf_aoa000 

fluent  2ddp  -t2  -pnmpi  -cnf=$PBS_NODEFILE  -g  -i 

journal_FLUENT_aoa000_152_tahoe . jou  >  FLUENT_aoa000_152 . out 

Some  output  is  directed  to  the  file  “out”,  while  any  errors  messages  that  occur  are 

directed  to  the  file  “error”.  While  only  1  node  is  used,  2  processors  on  that  node  are 

utilized.  To  utilize  only  1  processor,  simply  change  the  switch  “-t2”  to  “tl”  in  the  last 
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line  of  the  script.  The  directory  is  changed  to  that  of  where  the  PBS  script  is,  and  any 
previous  UDF  compilation  library  directory  is  removed.  Fluent®  is  started  in  2D  with 
double  precision  variables.  The  *.jou  journal  file  is  fed  to  Fluent®  to  be  executed  and  all 
screen  output  is  directed  to  a  *.out  file.  The  PBS  script  will  conclude  once  the  journal 
file  script  has  completed. 


# ! /bin/bash 

#PBS  -o  out 

#PBS  -e  error 

#PBS  -1  nodes=l : ppn=2 

#PBS  -j  oe 

Cd  $PBS_0_W0RKDIR 

numbers= '000 
+  02 
+  04 
+  06 
+  08 
+  10 
+  12 
+  14 
+  16  ' 

prenamein  = journal_FLUENT_aoa 
postnamein  =_152_tahoe . jou 
prenameout  =FLUENT_aoa 
postnameout=_l 52 . out 
libname  =  libudf_aoa 

for  i  in  $numbers 
do 

echo  $i  UDF  Library  being  Erased 
rm  -R  -f  $libname$i 

echo  $i  Started 

fluent  2ddp  -tl  -pnmpi  -cnf=$PBS_NODEFILE  -g  -i 

$prenamein$i$postnamein  >  $prenameout$i$postnameout 

echo  $i  Finished 
echo 

done 
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For  consecutive  serial  runs,  the  more  efficient  approach  was  performed.  This 
allows  for  the  simulations  to  be  called  one  at  a  time  for  the  entire  Angle  of  Attack  data 
set,  which  required  tremendously  less  amounts  of  user  interaction  between  runs  and 
hence  completes  the  set  faster.  However,  the  fastest  approach  is  to  run  all  of  the  set  in 
parallel,  assuming  the  processors  and  Fluent®  software  licenses  are  available. 

Journal  Script 

The  journal  file  is  the  heart  of  the  simulation’s  running.  The  only  input  into  this 
file  is  the  information  on  the  left  of  the  below  table.  To  help  with  understanding,  the 
prompts  that  appear  for  some  of  the  settings  are  shown  to  the  right. 

The  settings  for  Fluent®  reside  in  a  folder  system  that  can  be  accessed  via  text. 
To  move  up  one  folder  level  requires  the  command  “q”.  Moving  into  a  folder  can  be 
accomplished  by  typing  its  name  while  in  the  directory  it  exists  in. 

The  main  objective  of  the  journal  file  is  to  setup  the  simulation,  define  the  User 
Defined  Functions  (UDFs),  execute  the  simulation,  and  write  the  appropriate  data  for 
later  examination. 


First  we  read  in  the  Case  and  Data  File. 


file 

cd  file/ 

/file/ 

red 

Read  Case  File  &  Data 

NACAO  0  0  9_v0  7d6_base_FLUENT 

*.cas  and  *.dat  Case  File  &  Data 

respectively 

q 

cd  .. 

/ 
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The  user  defined  memory  is  setup  and  the  UDF  library  is  compiled. 


define 

cd  define/ 

/define/ 

user-defined 

cd  user-defined/ 

/define/user-defined/ 

user-def ined-memory 

Setting:  user- defined- memory 

5 

Number  of  User- Defined  Memory, 

UDM,  locations 

compiled- functions 

Setting:  compiled- functions 

compile 

load/unload/compile? 

Iibudf_aoa000 

Compiled  UDF  library  name 

yes 

Continue? 

temp_mom_src_trm_FLUENT . c 

Give  C- Source  file  names: 

First  file  name 

I!  II 

Next  file  name 

If  II 

Give  header  file  names: 

First  file  name 

compiled- functions 

Setting:  compiled- functions 

load 

load/unload/compile 

Iibudf_aoa000 

UDF  Library  Name 

Angle  of  Attack  set  at  0  degrees 

q 

cd  .. 

/define/ 
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Species  Transport  is  enabled  to  allow  for  the  fluid  Energy  Equation  to  be  used.  This  will 


allow  for  the  input  of  the  Temperature  UDF  Source,  but  is  not  required  for  the  input  of 
the  Momentum  UDF  Source. 


models 

cd  models/ 

/define/models/ 

viscous 

cd  viscous/ 

/define/models/viscous/ 

spalart-allmaras? 

Setting:  spalart-allmaras? 

yes 

Enable  the  Spalart-Allmaras 

Turbulence  model? 

q 

cd  .. 

/define/models/ 

species 

cd  species/ 

/define/models/species/ 

species-t ran sport  ? 

Setting:  species-transport? 

yes 

Enable  the  species  transport  model? 

mixture-template 

Select  an  available  mixture  material, 
(mixture-  template) 

q 

cd  .. 

/define/models/ 

q 

cd  .. 

/define/ 

The  boundary  conditions  for  the  velocity  inlet  are  set. 


boundary-conditions 

cd 

/define/boundary-conditions/ 

velocity-inlet 

Setting:  velocity- inlet 

velocity-inlet-6 

zone  id/name 

yes 

Velocity  Specification  Method: 
Magnitude  and  Direction? 

yes 

Reference  Frame: 

Absolute 

no 

Use  Profile  for  Velocity  Magnitude? 

15 . 2 

Velocity  Magnitude  (m/s) 
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no 

Use  Profile  for  X-Component  of  Flow 
Direction? 

+1 . 000000000000 

X-Component  of  Flow  Direction 
cos  (a) ,  a  =  Angle  of  Attack 

no 

Use  Profile  for  Y-Component  of  Flow 
Direction? 

+0 . 000000000000 

Y-Component  of  Flow  Direction 
sin  (a ) ,  a  =  Angle  of  Attack 

no 

Use  Profile  for  Temperature? 

288 . 15 

Temperature  (k) 

yes 

Turbulence  Specification  Method: 
Modified  Turbulent  Viscosity 

no 

Use  Profile  for  Modified  Turbulent 
Viscosity? 

0 . 001 

Modified  Turbulent  Viscosity  (m2/s) 

no 

Use  Profile  for  h2o  mass  fraction? 

0 

h2o  mass  fraction 

no 

Use  Profile  for  o2  mass  fraction? 

0 

o2  mass  fraction 

The  boundary  conditions  for  the  fluid  are  set.  The  UDF  Source  Terms  for  Energy,  X  & 


Y  Momentum  are  tied  in  as  well. 


fluid 

Setting:  fluid 

fluid 

zone  id/name 

yes 

Specify  source  terms? 

no 

Use  Constant  Mass  (kg/m3-s)  source? 

no 

Use  UDF  for  Mass  (kg/m3-s)  source? 

no 

Use  Constant  X  Momentum  (n/m3) 
source? 

yes 

Use  UDF  X  Momentum  (n/m3)  source? 

"x_momentum_source : : libudf_aoa000 " 

udf-  name 
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no 

Use  Constant  Y  Momentum  (n/m3) 
source? 

yes 

Use  UDF  Y  Momentum  (n/m3)  source? 

"y_momentum_source : : Iibudf_aoa000 " 

udf-  name 

no 

Use  Constant  Modified  Turbulent 
Viscosity  (kg/s2-m)  source? 

no 

Use  UDF  Modified  Turbulent  Viscosity 
(kg/s2-m)  source? 

no 

Use  Constant  h2o  (kg/m3-s)  source? 

no 

Use  UDF  h2o  (kg/m3-s)  source? 

no 

Use  Constant  o2  (kg/m3-s)  source? 

no 

Use  UDF  o2  (kg/m3-s)  source? 

no 

Use  Constant  Energy  (w/m3)  source? 

yes 

Use  UDF  Energy  (w/m3)  source? 

"temp_source : : Iibudf_aoa000 " 

udf-  name 

no 

Specify  fixed  values? 

yes 

Motion  Type:  Stationary? 

0 

X- Origin  of  Rotation- Axis  (m) 

0 

Y- Origin  of  Rotation- Axis  (m) 

no 

Deactivated  Thread 

no 

Laminar  zone? 

no 

Porous  zone? 

q 

cd  .. 

/define/ 

Setting  the  material  properties  of  the  fluid,  air. 


materials 

cd  materials/ 

/define/materials/ 

change -ere ate 

Setting:  change-create 

air 

material-  name> 
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air 

material  name 

yes 

air  is  fluid 
change  Density? 

constant 

Density  methods: 
new  method 

0 . 993 

value  (kg/m3) 

no 

change  Cp  (Specific  Heat)? 

no 

change  Thermal  Conductivity? 

yes 

change  Viscosity? 

constant 

Viscosity  methods: 
new  method 

1 . 70e-5 

value  (kg/m-s) 

no 

change  Molecular  Weight? 

no 

change  L-J  Characteristic  Length? 

no 

change  L-J  Energy  Parameter? 

no 

change  Thermal  Expansion 

Coefficient? 

no 

change  Degrees  of  Freedom? 

no 

change  Speed  of  Sound? 

q 

cd  .. 

/define/ 

Setting  the  operating  pressure  of  the  simulation. 


operating-conditions 

cd  operating- conditions/ 
/define/operating-  condition  s/ 

operating-pressure 

Setting:  operating- pres  sure 

78669 

operating  pressure  (pascal) 

q 

cd  .. 

/define/ 

q 

cd  .. 

/ 
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To  increase  the  speed  at  which  the  simulations  could  be  performed,  parallelization  on  2 
processors  was  implemented.  This  snippet  partitions  the  grid. 


parallel 

cd  parallel/ 

/parallel/ 

partition 

cd  partition/ 

/parallel/partition/ 

auto 

cd  auto/ 

/parallel/partition/auto/ 

use-case- file-met hod 

Setting:  use-case- file- method 

yes 

use  case- file  partition  method? 

q 

cd  .. 

/parallel/partition/ 

q 

cd  .. 

/parallel/ 

q 

cd  .. 

/ 

To  ensure  accurate  reporting  results  for  the  coefficients  of  lift  and  drag,  reference  values 
need  to  be  set  correctly. 


report 

cd  report/ 

/report/ 

reference-values 

cd  reference- values/ 
/report/reference-values/ 

compute 

cd  compute/ 

/report/reference-values/compute/ 

velocity-inlet 

Setting:  velocity- inlet 

velocity-inlet-6 

zone  id/name 

q 

cd  .. 

/report/reference-values/ 

area 

Setting:  area 

0 .202 

reference  area  (m2) 

q 

cd  .. 

/report/ 

q 

cd  .. 

/ 
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The  correct  time  step  needed  to  be  set. 


solve 

cd  solve/ 

/solve/ 

set 

cd  set/ 

/solve/set/ 

time-step 

Setting:  time- step 

0 . 001 

time  step  (s) 

q 

cd  .. 

/solve/ 

Output  of  the  coefficient  of  lift  to  a  file  is  setup. 


monitors 

cd  monitors/ 

/solve/monitors/ 

force 

cd  force/ 

/solve/monitors/force/ 

lift-coefficient 

Setting:  lift- coefficient 

yes 

monitor  cl? 

4 

zone  id/name(l) 

0 

zone  id/name(2) 

no 

print  cl  data? 

yes 

write  cl  data? 

"cl-history_aoa000 " 

cl  data  file  name? 

no 

plot  cl  data? 

no 

plot  per  zone? 

+0 . 000000000000 

x- component  of  lift  vector 
sin  (a ) ,  a  =  Angle  of  Attack 

+1 . 000000000000 

y- component  of  lift  vector 
cos  (a) ,  a  =  Angle  of  Attack 
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Output  of  the  coefficient  of  drag  to  a  file  is  setup. 


Setting:  drag- coefficient 
monitor  cd? 


drag-coefficient 

Setting:  drag- coefficient 

yes 

monitor  cd? 

4 

zone  id/name(l) 

0 

zone  id/name(2) 

no 

print  cd  data? 

yes 

write  cd  data? 

"cd-history_aoaOOO " 

cd  data  file  name? 

no 

plot  cd  data? 

no 

plot  per  zone? 

+1 . 000000000000 

x- component  of  lift  vector 
cos  (a) ,  a  =  Angle  of  Attack 

+0 . 000000000000 

y- component  of  lift  vector 
sin  (a ) ,  a  =  Angle  of  Attack 

q 

cd  .. 

/solve/monitors/ 

q 

cd  .. 

/solve/ 

initialize 

cd  initialize/ 

/solve/initialize/ 

Default  reference  values  for  pressure  and  velocity  were  setup  and  the  flow  initialized 

with  these  values. 

set-defaults 

cd  set- defaults/ 

/solve/initialize/  set- defaults/ 

pressure 

Setting:  pressure 

0 

Default  value  for  Gauge  Pressure 

x-velocity 

Setting:  x-velocity 

0 

Default  value  for  X  Velocity 

y-velocity 

Setting:  y-velocity 
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0 

Default  value  for  Y  Velocity 

q 

cd  .. 

/solve/initialize/ 

initialize-flow 

Setting:  initialize-flow 

q 

cd  .. 

/solve/ 

q 

cd  .. 

/ 

Now,  just  prior  to  starting  the  iterative  time  solving,  the  weighting  function  is  called. 


define 

cd  define/ 

/define/ 

user-defined 

cd  user- defined/ 

/define/u  ser-  defined/ 

execute -on-demand 

Setting:  execute-on-demand 

" cell_weight_on_demand : : libudf_aoa000 " 

Execute  on  demand  function  name 

q 

cd  .. 

/define/ 

q 

cd  .. 

/ 

Start  solver  for  1500  time  steps  with  a  maximum  of  20  iterations  per  time  step. 


solve 

cd  solve/ 

/solve/ 

dual-time- iter ate 

Setting:  dual- time- iterate 

1000 

Number  of  physical  time  steps 

20 

Number  of  iterations  per  time  step 

q 

cd  .. 

/ 
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Once  solution  has  completed  all  of  the  time  steps,  write  the  Case  and  Data  to  a  file  for 


later  review. 


file 

cd  file/ 

/file/ 

wed 

Write  Case  &  Data  File 

NACAO  0  0  9_v0  7d6_aoa0  0  0_FLUENT_end 

case/data  file  name 

yes 

OK  to  overwrite? 

q 

cd  .. 

/ 

Exit  the  program. 


exit  program 
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